Introduction
Gene content is the number of copies of a particular allele in a genotype of an animal (0, 1 or 2). Information on gene content can be used to estimate the effect of a candidate gene, in evaluation of the total genetic merit including polygenic and major gene effects or for identifying favourable selection candidates based on disease gene status. Gene content is known for genotyped animals but usually only a fraction of animals has their genotype known. For ungenotyped animals gene content can be predicted from the records on their typed relatives. Methods for exact calculation of gene content are not practical for large pedigrees (Lauritzen and Sheehan, Reference Lauritzen and Sheehan2003). Approximate methods have also been proposed but they are complex, require custom computer programs and may be slow when applied to pedigrees of millions of animals (Kong, Reference Kong, Keramidas and Kaufman1991; Thallman et al., Reference Thallman, Bennett, Keele and Kappes2001).
The double-muscling condition is known in several breeds. The first description was made in Shorthorn cattle and it seems that the Belgian Blue breed (BBB) has inherited the condition from them. After 1950 the selection of double-muscled animals in the BBB was increased strongly and selection for milk production discontinued. In 1974 two strains were defined: meat and dual-purpose. The milking cows are mostly dual-purpose registered (DP-BBB), however some meat strain registered animals (M-BBB) are still milked. The double-muscling condition is under the control of a major gene but modified by other genes. Charlier et al. (Reference Charlier, Coppieters, Farnir, Grobet, Leroy, Michaux, Mni, Schwers, Vanmanshoven, Hanset and Georges1995) identified the major gene in BBB as being a mutation inactivating myostatin. The M-BBB animals are homozygous for the mutated allele called muscular hypertrophy (mh). In DP-BBB three genotypes are encountered: +/+, mh/+ and mh/mh. All DP-BBB sires have to be genotyped before they can be registered and over the years top breeders genotyped their cows too. Breeders use the genotype at the myostatin locus as a criterion in selection.
In this paper we present a practical method to calculate gene content in large pedigrees. In simulation study we compare the proposed method with other approximate techniques: Markov chain Monte Carlo (MCMC), iterative peeling and the procedure of Israel and Weller (Reference Israel and Weller1998). We consider the calculated gene content for further calculation of approximate genotype probabilities. Finally, we exemplify the use of the method in a study of the mh mutation for the myostatin locus in DP-BBB cows under milk recording.
Material and methods
Method
To present the method we consider inheritance at single biallelic locus with hypothetical alleles B and b. Let q be the number of copies of allele B carried by an individual (gene content). For an animal that has been genotyped q is 0, 1 or 2. For an individual that has not been genotyped but is the progeny of genotyped parents the expected value for q is E(q p) = 0.5(q s+q d), with q s, q d and q p being the content of allele B for sire, dam and progeny. We can replace q by its deviation from population mean, d = q − μ. Then, expected value for a progeny of genotyped parents is E(q p) = μ+0.5(d s+d d). If mother has not been genotyped her expected value has to be replaced by population mean, q d = μ. Further, if a relative with genotype known is not a parent or offspring expected vales for q can be derived using a recursive approach. For example, if maternal grand-sire (mgs) has been genotyped, the expected value of q for mother is E(q d|d mgs) = μ+0.5d mgs, and expectation for progeny is E(q s|d s, d mgs) = μ+0.5d s+0.25d mgs. Up to this point the method is similar to the approach by Israel and Weller (Reference Israel and Weller1998). Their method, however, ignores information on descendents that usually are the main source of information. To use information from all genotyped relatives we apply known techniques developed for continuous traits. Further, we treat gene content as a continuous variable and assume that the relationship between gene contents is linear, and that the covariance between gene contents is proportional to the additive relationship between animals. Although the assumptions are violated we will show in simulation study that the method is robust enough to enable practical calculations. Under these assumptions the conditional expectation of gene contents for all ungenotyped animals given molecular and pedigree data is
where qx is an unknown vector of gene contents for animals that have not been genotyped, 1 is a vector of one's, qy is a known vector of gene contents for animals that have been genotyped, Axy is an additive relationship matrix between individuals with unknown genotype and their genotyped relatives, Ay is the additive relationship matrix for genotyped individuals. Note that the mean value (μ) can be easily calculated if we assume that animals are typed without error.
A more practical way of computation is derived under an incomplete penetrance model. Here, the incomplete penetrance model incorporates the probability of errors in marker phenotypes and is more accurate for real marker data. Assuming small error variance component () in the variability of q, we use mixed model equations that provides at the same time the solution for μ:
where dy is a vector of gene content deviations for animals with genotype records, dx is a vector of gene content deviations for unobserved animals, M is incidence matrix linking qy to that can be rewritten as A is the additive relationship matrix of the structure The order of animals was chosen to be genotyped before ungenotyped, in practice however, the order is of user choice. For this model the assumptions of BLUP methodology are applied. Some small error variance is required to solve the system of equations. In practice, for the model of almost complete penetrance (small probability of mistyping) the ɛ can be set equal to 0.01 or smaller. A large value for ɛ would mean that probability of mistyping or pedigree error is high. Such data would be useless.
This method has some important advantages: standard genetic evaluation software can be used and gene contents can be easily computed for the whole population; population mean gene content will be computed at the same time; the model can be adapted to account for genetic groups due to different origins of animals; ratio ɛ between variances can be small but still different from zero, therefore the method allows for rare but not impossible genotyping errors, ɛ may also account for pedigree errors, that could be in certain cases up to over 20% of the animals (Banos et al., Reference Banos, Wiggans and Powell2001).
Simulation study
The proposed method has been compared with three other methods: MCMC approach, iterative peeling and the method developed by Israel and Weller (Reference Israel and Weller1998). All three methods approximate genotype probabilities for each animal in a pedigree. Exact genotype probabilities are marginal probabilities obtained by summing over all possible genotype configurations which are consistent with observations, weighting the configurations by their probability of occurrence as calculated under the assumed genetic model. On a small pedigree exact genotype probabilities can be calculated using one of the graphical models and reverse or simultaneous peeling (Thompson, Reference Thompson1981). The problem considered in the simulation study was too complex for exact calculation, therefore only approximate methods were applied. The approximate probabilities computed by MCMC, iterative peeling and the method of Israel and Weller (1998) were used to calculate expected gene content of an animal applying the following formula E(q) = 2P(B B)+P(B b). These results were used to compare with values obtained from the proposed method. The proposed method was implemented with BLUPF90 (Misztal, Reference Misztal1999) and all other methods with our own software. To compare the methods 200 data sets have been randomly generated for a real bovine pedigree and the data sets have been analysed by all four methods.
The pedigree used in the simulations was extracted from the DP-BBB pedigree data. All cows of a single herd and their ancestors were represented. The pedigree consisted of 1082 individuals: 190 founders, 892 non-founders. Within the pedigree 716 nuclear families (pairs of parents with their siblings) and 531 loops were identified. On a pedigree graph, a loop occurs when an animal can be connected to itself through parents or progeny.
The milking BBB population is partially genotyped at the myostatin locus (recent actual DP-BBB animals). The genotype records were available for 129 individuals born from 1969 through 2001 (Table 1). The structure of this real data set was used in the simulation but real records were replaced by simulated values. Simulation of a complete data set for all pedigree members enables a comparison of the calculated gene content for an animal to its simulated value. A complete data set was generated for biallelic polymorphic system by the use of random process of gene dropping. To make the structure of the simulated data set similar to real data set, the complete simulated data set was reduced to 129 records for the animals, for which real records exist (Table 1), and the reduced data set was analysed by all four methods. Two frequencies of B allele were considered: 0.4 and 0.2. One hundred data sets were generated and analysed under each of the two gene frequencies.
The accuracy of the proposed method for different proportions of missing genotypes was tested on huge bovine pedigree of 907 903 animals (see below). Genotypes at a single biallelic locus were simulated under equal allelic frequency by the gene dropping method. Three proportions of missing genotypes were considered: 90%, 70% and 50%.
MCMC method
In the MCMC method genotype probabilities are calculated from multiple random samples of animals' genotypes which are consistent with pedigree and molecular information. To sample the genotype configuration over the pedigree we used the whole locus sampler (Kong, Reference Kong, Keramidas and Kaufman1991). The sampler uses a peeling algorithm to sample genotypes of all animals jointly from the desired distribution. The pedigree was too complex to peel the pedigree considering the total genotype space. Thus, the genotype space was reduced using the concept of set-recoding and fuzzy inheritance (O'Connell and Weeks, Reference O'Connell and Weeks1995). In brief, two possible alleles of an individual were replaced with a set s = {B,b} if they were not observed on the individual nor on his descendents. For example, if within a family of sire, dam and progeny only the sire has its genotype known, say BB, their recoded ordered genotypes were {B}|{B}, {B,b}|{B,b} and {B}|{B,b}, respectively. Considering fuzzy inheritance, if a parent has genotype s 1|s 2 and its child has allele set s 3, then s 3 is inherited from the parent if s 1 ⊆ s 3 or s 2 ⊆ s 3. To sample genotypes we used the following scheme in each iteration: (1) a random ordered genotype configuration on reduced genotype space was obtained by random propagation after peeling the pedigree toward a root, the ordered genotypes were considered to differentiate Bb and bB heterozygotes, (2) segregation indicators (SI) were set based on current ordered genotypes or sampled randomly if both values were possible, (3) the gene frequencies were sampled based on current allelic types (not recoded) in founders, (4) each recoded allele of a founder was replaced by random allele based on current allelic frequencies, and, starting from the top of the pedigree, each recoded allele of a non-founder was replaced by parental allele according to its SI value. Note, for these parts of the pedigree, for which no data was available, the procedure was similar to simple gene dropping. This sampling procedure guarantees irreducibility and has good mixing properties. When experimenting with small pedigrees it provided solutions almost equal to exact values calculated by the use of a peeling algorithm (Elston and Stewart, Reference Elston and Stewart1971) and PAP package (Hasstedt, Reference Hasstedt2002). For each of the 200 simulated data sets we ran 100 000 iterations collecting samples from all iterations. The final solutions for genotype probabilities were calculated as averages from 100 000 samples. Note, as the whole pedigree was sampled jointly the samples were correlated only through gene frequency and the autocorrelation of the samples was very low.
Iterative peeling
The iterative peeling exploits the fact that given genotype probabilities of an animal's neighbours (parents, offspring and mates) the genotype of that animal depends only on neighbours and it is independent from the rest of the pedigree. Starting from some initial values genotype probabilities are updated individual by individual until convergence (Van Arendonk et al., Reference Van Arendonk, Smith and Kennedy1989; Janss et al., Reference Janss, Van Arendonk and Van der Werf1995; Wang et al., Reference Wang, Fernando, Stricker and Elston1996, Thallman et al., Reference Thallman, Bennett, Keele and Kappes2001). It was demonstrated on small pedigrees that the iterative peeling may lead to genotype probabilities close to exact calculations (Fernandez et al., Reference Fernandez, Fernando, Guldbrandtsen, Totir and Carriquiry2001), however, the behaviour of iterative peeling on large complex pedigrees is unknown. We used the method for comparison with the proposed approach on simulated and real data sets. In simulation study the gene frequency of mh was assumed to be known at the simulated value, and for analysis of real data on the myostatin locus the frequency was set equal to 0.05.
Analysis of myostatin locus
Sparse molecular data on the myostatin locus was used to exemplify the utilisation of the method on a huge bovine pedigree for which exact method is impossible and MCMC methods difficult. The population consisted of 907 903 animals born between 1960 and 2005 derived from Belgian dairy and dual-purpose cattle data base. Within the population 34% individuals were Holstein Friesian, 26% were BBB under milk recording including both M-BBB and DP-BBB strains, 10% were Red and White, 1% of other origin and 29% were crosses of various breeds. Within the population 1865 animals have been genotyped for the mh mutation at the myostatin locus. Additionally, recent M-BBB sires were considered homozygous for the mh gene and all dairy breeds were assumed free of the mh gene. In total, there were 381 742 records on the mh gene content and the mh gene content was estimated for the rest (58%) of the population. Including dairy breeds was necessary because ancestors could come from those breeds. The data set was also analysed with modified model linking progeny of unknown parents to unknown-parent groups (Westell et al., Reference Westell, Quaas and Van Vleck1988).
Results
Simulation study
The study showed that calculation of gene content when only 12% of animals have genotype records may be difficult. Correlation coefficients between simulated and calculated gene contents under gene frequency 0.4 and 0.2 were 0.50 and 0.47, respectively. Similar correlations were calculated when gene contents were calculated by the use of MCMC approach (0.52 and 0.42) and iterative peeling (0.52 and 0.40). The proposed method showed significant improvement over the method of Israel and Weller (Reference Israel and Weller1998), for which the correlation coefficients were 0.38 and 0.38.
Considering MCMC solutions and allelic frequency of 0.4, the mean of calculated gene content for animals with simulated genotype bb (simulated q = 0) was 0.59, for alternative homozygote BB (simulated q = 2) only 1.14. The range between these means was only 28% of the true range of 2 (Table 2). For lower degree of polymorphism (allelic frequency of 0.2) the range between means was 26% of the true range. In general, the bias increases for rare genotype. Standard deviation for rare genotype BB was smaller than for more common genotypes, a consequence of the strong bias of the mean for this genotype. Iterative peeling provided solutions close to MCMC values. The proposed method behaves slightly worse. For population allelic frequency of 0.4, the range between means for alternative homozygotes was 25% of the true range, and for frequency of 0.2 the range decreased to 21%. In the later case the bias decreased for the common homozygote and increased for the rare homozygote. We considered the errors defined as the difference between simulated and calculated gene contents. Compared with the MCMC method the mean square error (results not shown) was smaller in case of the more common genotype (bb) and higher for the rare genotype (BB). For allelic frequency of 0.2 the total mean square error calculated across three genotypes was considerably lower (.222) than for MCMC method (.288).
The behaviour of the proposed method for different proportions of missing genotypes is presented in Table 3. For comparison, the data sets were also analysed by the use of the iterative peeling method. For the proposed method, the correlation between simulated and estimated gene content was 0.54 when 90% of animals had no genotype records, and increased to 0.62 and, 0.67 for the data sets with 70 and 50% of missing genotypes in the pedigree. The corresponding coefficients calculated for the iterative peeling method were slightly lower: 0.48, 0.58 and 0.65. The differences between means for alternative homozygotes were comparable between the two methods, only for the most dense data set (50% of missing genotypes) did the iterative peeling algorithm slightly surpassed the proposed method. The amount of time required to complete the calculation by the use of the proposed method was much shorter than that needed for the iterative peeling algorithm. In the case of 90% missing genotypes the proposed method took 4 min (by the preconditioned conjugate gradient algorithm) and the iterative peeling method took almost 4 h (32 iterations). For the proposed method the amount of time increased with the number of records but rapidly decreased in the case of the iterative peeling. For the dense data set (50% of missing genotypes) the proposed method was still three times faster than the iterative peeling algorithm.
Analysis of myostatin locus
Content of the mh gene was estimated for 526 161 animals, including 235 133 animals of DP-BBB breed. Two models were used for the proposed method: the standard model as presented above and a model with unknown parent groups as described by Westell et al. (Reference Westell, Quaas and Van Vleck1988). The model with unknown parent groups was used to account for selection for mh allele in BBB animals. Cows were assigned to six groups: five for BBB born in different time periods and one for all non-BBB animals. Similar groups were created for sires. For DP-BBB animals mean mh content was calculated for each year of birth (Figure 1). Results from the two models showed rapid increase in the frequency of mh allele. The standard model indicates that the frequency of mh gene increased since 1970, while the model with unknown parent groups shows that a significant increase of the gene started 10 years later. It is unclear which of the two models better fits the history of the DP-BBB breed.
For comparison with the proposed method, the content of mh gene was also calculated by the use of the iterative peeling algorithm (Figure 1). If only one genetic group was assumed, the calculated mean mh content for old animals was very high. This result should be considered incorrect because it is known that in the early years the frequency of the mh gene was rather low. An additional analysis was performed by the use of the iterative peeling and multiple genetic groups. All founders were assigned to the same 12 groups and frequency of mh gene in each group was assumed to be known. The assumed values for the mh frequency in each group were taken from the results obtained by the proposed method. Two analyses by the use of the iterative peeling method under multiple genetic groups model were performed, one with the mh gene frequency based on the results from the proposed method under standard model, and the other under the Westell et al. (Reference Westell, Quaas and Van Vleck1988) model. The results of these two analyses were similar and little influenced by the assumed mh frequency. The results from the iterative peeling were similar to the results calculated by the proposed method under the standard model, even if the mh frequency was based on the Westell et al. (Reference Westell, Quaas and Van Vleck1988) model.
Discussion
There are two approaches to approximate gene content for the ungenotyped animals in huge pedigrees using the information from all their genotyped relatives: the Monte Carlo approach and the iterative peeling approach. The first is more flexible and therefore can be more useful in practice. The problem with these methods is that it is difficult to prove that they provide good solutions in large pedigrees unless the pedigree is peelable and exact values can be calculated. There are number of possible Monte Carlo samplers, however, all of them are Markov chain samplers unless again the pedigree is peelable. Advanced samplers use the peeling or the iterative peeling computation to sample from desired probability distribution (Kong, Reference Kong, Keramidas and Kaufman1991; Fernandez et al., Reference Fernandez, Fernando, Guldbrandtsen, Totir and Carriquiry2001). MCMC approaches are very useful in the analysis of genetic data on pedigrees, however, the construction of efficient Markov chain and monitoring its convergence requires extensive additional experience from the user of these methods. Moreover, both advanced MCMC samplers and the iterative peeling require relatively sophisticated computer programs and availability of such algorithms is limited. In this paper we propose a practical method to approximate the gene content in large animal pedigrees. All computations can be done using one of a few trusted statistical packages for genetic evaluation.
The time needed to complete computation depends on the quality of the software and the convergence criteria. In case of the MCMC method the time also depends on the mixing behaviour and the required precision of estimation. In practice, limited amount of time can be allocated to a method to do computation. Under time limit, some MCMC methods may provide inaccurate solutions (Totir et al., Reference Totir, Fernando, Dekkers, Fernandez and Guldbrandtsen2003). The proposed method uses a simplified genetic model and the simplification drastically reduces the amount of time required to complete calculation. In contrast to the alternative methods, the amount of time does not increase with the number of missing genotypes. The capability of the model to calculate unknown gene content is comparable with the models usually used by the MCMC and iterative peeling approaches. The proposed method can be considered practical for large animal pedigrees with a high proportion of missing genotypes.
Some genetic problems would require genotype probabilities. In contrast to other methods, the proposed method does not provide genotype probabilities. Approximate genotype probabilities can be derived from the calculated gene content. For this we treat half the gene content as probability of B allele and calculate the genotype probability from Hardy-Weinberg proportions. Note that the crucial assumption, which should hold for this calculation is that of equal viability of all gametes and zygotes, and this assumption has already been made when deriving the expected value for q from the mean of parental values. We believe that this assumption will hold in most cases unless a gene causes death at early age. An additional assumption required is that of random mating in base population. We compared genotype probabilities calculated in this way to genotype probabilities obtained directly from the iterative peeling method applied to a model with multiple genetic groups. For this comparison a population of DP-BBB cattle (235 133 animals) was considered. Note, within the population only 1284 animals had their genotype records, however, many DP-BBB animals have ancestors in dairy breeds, which were assumed free of the mh allele, and some were related to M-BBB animals which were considered all mh/mh. This fact was also taken into account during computations. The results of the comparison are presented in Table 4. The standard deviation of the absolute difference between alternative calculations was less than 0.12, however, for a few animals large dissimilarity between results have occurred. The highest difference between two alternative solutions was 0.95. The high dissimilarity occurred usually for an animal that had no genotype record but its only possible genotype is fixed based on typed relatives.
The proposed approach does not detect inconsistency between genotype and pedigree records. However efficient algorithms to detect non-mendelian records exist and they can be applied before calculation. One of possible options is to use a genotype elimination algorithm as proposed by Lange and Goradia (Reference Lange and Goradia1987). The algorithm works iteratively on the nuclear families and has enough power to detect most of the errors that would result in null likelihood under a complete penetrance model. Application of the genotype elimination algorithm may add some data by concluding the only possible genotype for an animal.
The proposed method is similar to the best linear unbiased prediction applied to discrete data. Application of BLUP methodology to discrete variables has been considered in the context of breeding value estimation (Gianola, Reference Gianola1980; Hoeschele, Reference Hoeschele1986). The main problem in the analysis of a discrete trait is that the restriction that the sum of response probabilities across all categories must total 1 is not taken into consideration. Numerous generalised linear mixed models have been proposed to study underlying continues variables from discrete data. These models may improve precision of the proposed method. Unlike the mixed model equations, however, the estimation equations for generalised linear mixed model must be solved iteratively. Such a method would not be more practical than the iterative peeling method. Moreover, the assumption of the linear relationship between gene contents would remain violated.
In the analysis of mh gene we considered multiple genetic groups to demonstrate that a genetic model can be easily adapted for a particular pedigree. The proposed method uses the general methodology of the linear mixed models and the trusted computer applications covering a wide range of various linear models are easily available. Adaptations of the model, as separate means of gene content (μ) or, as shown in this paper Westell groups can be easily specified for different breeds or birth years of animals (Westell et al., Reference Westell, Quaas and Van Vleck1988). Although, a model used by the alternative approaches can also be adapted to fit particular data, the available computer programs for these algorithms are usually hard to modify.
The proposed method does not require that the frequency of alleles among founders are known. In the analysis of mh gene it was demonstrated that this assumption is important for the iterative peeling approach. Certainly, the effect of incorrect assumption on gene frequency increases with the proportion of the missing genotypes. Again, this indicates that the proposed method will be especially practical for large pedigree with sparse data.
The proposed method has a property that may be considered useful when the analysed data contain hidden errors. The example is a sire with hundred daughters, all BB but one bb, which suggests a typing or pedigree error but yields positive likelihood under complete penetrance model. The exact genotype probability for the sire is P(Bb) = 1 and consequently q = 1. The presented algorithm is robust against this potential error in the sense that it calculates a value depending on the number of BB progeny, here q = 1.98. This property is useful in large animal pedigrees because in most cases they are not subjected to close scrutiny. Accounting for errors in the iterative peeling approach would require relaxation of the assumed genetic model.
In future, genetic evaluation will more often be based on the identified polymorphisms affecting directly quantitative traits of interest. It can be anticipated that usually biallelic genes as single nucleotide polymorphisms (SNP) will be identified. The proposed method is well fitted to this kind of polymorphism. In case of a multiallelic polymorphism, however, the application of the proposed method is still possible by estimating the content of each allele separately. The method cannot be easily modified to use information from linked markers.
It was also shown that sparse data provides little information on true gene content of an ungenotyped animal, and consequently, such data would not contribute much to identify favorable selection candidates. Also estimation of candidate gene effects will be difficult. Accurate prediction of gene content would require typing more individuals. The way the animals are chosen for genotyping may influence the estimate of candidate gene effects (Israel and Weller, Reference Israel and Weller1998). Phenotypes under the influence of the typed gene could contribute additional information on gene content.
The proposed method can also be useful in initialisation of MCMC genotype sampling for solving some more complex problem, e.g. a multilocus problem. Single site and blocked samplers require initial genotype configuration consistent with observed data. A consistent starting configuration may be obtained by ‘gene-dropping’, drawing the genes for the founders from some assumed base population and sampling subsequent generation according to Mendelian segregation law, rejecting configurations which are inconsistent with data. For large pedigrees, however, gene-dropping has low acceptance ratio. The presented method can improve the acceptance ratio of gene-dropping algorithm by including approximate gene contents in calculation for each sampled genotype.
In this paper we presented a practical method to calculate gene content for ungenotyped animals of large and complex pedigrees. The method is development out of the approach by Israel and Weller (Reference Israel and Weller1998) and can be used as alternative to more advanced approaches. The method can be easily modified to correctly use different base populations. Approximate genotype probabilities can also be calculated. The calculated values can be used to estimate the additive effect of a candidate gene or to support decision in marker-assisted selection.
Acknowledgements
Nicolas Gengler, who is Research Associate of the National Fund for Scientific Research (Brussels, Belgium), acknowledges this support. Additional support was provided through grant 2.4507.02F (2) of the National Fund for Scientific Research. Manuscript review by Dr George R. Wiggans (Animal Improvement Programs Laboratory – ARS – USDA, Beltsville; USA) and the support of the Walloon Breeding Association (AWE) and the Walloon Regional Ministry of Agriculture (project D31-1112) are gratefully acknowledged.