Inbreeding and Purging
Several measures of inbreeding and purging can be computed, based on the probability of allele identity by descent of individuals of the pedigree. All functions related to inbreeding and purging are prefixed with ip_.
Wright’s inbreeding coefficient
The inbreeding coefficient (\(F\), Wright 1922), here also referred to as standard inbreeding, is defined as the probability that an individual inherits two alleles derived from the same ancestor (i.e. identical by descent, IBD). In pedigreed populations, this can be calculated for an individual \(i\) as the kinship coefficient of its parents \(j\) and \(k\) (\(F_{i} = f_{j,k}\)), which can be calculated as:
\[f_{j,k {} (j=k)}=\frac{1}{2}(1+F_{j})\]
\[f_{j,k { } (j\neq k)}=\frac{1}{2}(f_{j,k_{d}}+f_{j,k_{s}})\] Where \(k_{d}\) and \(k_{s}\) refer to \(k\)’s dam and sire (see Falconer & Mackay 1996).
The function ip_F computes the inbreeding coefficient, given an input pedigree. Note that the value of \(F\) will be saved in a new column of the dataframe, as it is usually convenient to save it this way to simplify the computation of further inbreeding and purging parameters, as well as for later analyses.
The example below shows the inbreeding coefficient of William E. Darwin (son of Charles R. Darwin and Emma Wedgwood).
darwin <- darwin %>% purgeR::ip_F()
darwin %>% dplyr::filter(names == "William Erasmus Darwin")
#> id dam sire names Fi
#> 1 60 54 52 William Erasmus Darwin 0.06298828\(F\) can also be estimated based on population estimates of the effective population size \(N_{e}\) and generation numbers, using the classical expression (Falconer and Mackay 1996):
\[F_{t} = 1 - (1-\frac{1}{2N})^{t}\] This can be achieved with the function exp_F (e.g. exp_F (Ne = 50, t = 50)).
Partial inbreeding coefficient
As mentioned above, IBD happens when alleles are inherited from the same ancestor and appear in homozygosis. Thus, \(F_{i}\) can be partitioned as the additive contribution of its ancestors to \(F_{i}\). The partial inbreeding coefficient \(F_{i(j)}\) is defined as \(i\)’s probability of IBD for alleles coming from ancestor \(j\). It can be computed from partial kinship coefficients (\(f_{p_{1},p_{2}(j)}\), where \(p_{1}\) and \(p_{2}\) refer to \(i\)’s parents), so that \(F_{i(j)}=f_{p_{1},p_{2}(j)}\), using the tabular method as described by Gulisija & Crow (2007). Given an ancestor \(j\):
- All \(f_{p_{1},p_{2}(j)}\) values are initialized to \(0\), except for the diagonal entry corresponding to founder \(j\), which takes value \(1/2\).
- Values for intermediate ancestors are computed as follows:
- When \(p_{1}=p_{2}\): \(f_{p_{1},p_{2}(j)} = f_{p_{1},j}+\frac{1}{2}f_{p_{1d},p_{1s}}\), where \(p_{1d}\) and \(p_{1s}\) are \(p_{1}\)’s parents.
- When \(p_{1}\neq p_{2}\): \(f_{p_{1},p_{2}(j)} = \frac{1}{2}(f_{p_{1d},p_{2}}+f_{p_{1s},p_{2}})\), where \(p_{2}\) is older than \(p_{1}\).
The function ip_Fij will return a matrix object with all possible values of the partial inbreeding coefficient. In that matrix, the value in row \(i\) and column \(j\) indicates the probability of IBD of individual \(i\) for alleles coming from ancestor \(j\). Values in the upper diagonal of the matrix always take values of zero. Of course, the summation of \(F_{i(j)}\) over every column \(j\) equals \(F_{i}\) when \(j\) are founder ancestors.
m <- ip_Fij(arrui, mode = "founders") # ancestors considered are founders (by default)
base::rowSums(m) # this equals ip_F(arrui) %>% .$FiBy default, ip_Fij only considers partial inbreeding conditional to founders, but it can also be extended to any ancestor using the mode = "all" argument. A custom number of individuals can also be used (see ?ip_Fij). Note however that for a large number of individuals, the computation of this matrix may require a substantial amount of time. In every case, columns of the returned matrix are sorted by ancestor identity use.
Figure below shows the contribution of the two founders in the Barbary sheep pedigree to inbreeding values \(F>0.35\).
arrui <- arrui %>% purgeR::ip_F()
tibble::tibble(founder1 = m[, 1], founder2 = m[, 2], Fi = plyr::round_any(arrui$Fi, 0.025)) %>%
tidyr::pivot_longer(cols = c(founder1, founder2), names_to = "Founder", values_to = "Fij") %>%
dplyr::group_by(Fi, Founder) %>%
dplyr::summarise(Fij = sum(Fij)) %>%
ggplot() +
geom_bar(aes(x = Fi, y = Fij, fill = Founder), stat = "identity", position = "fill") +
scale_x_continuous("Inbreeding coefficient (F)", limits = c(0.35, 0.625)) +
scale_y_continuous("Partial contribution to F (in %)", labels = scales::percent_format()) +
scale_fill_manual(values = c("darkgrey", "black")) +
theme(
panel.background = element_blank(),
legend.position = "bottom"
)
Alternatively, partial inbreeding can also be estimated via genedrop simulation (using option genedrop). This will however result in less precise estimation of \(F_{i(j)}\), and might only be convenient to use in terms of performance for very large and complex pedigrees. In these cases, a value of genedrop = 100 might give results that are well correlated with exact estimates (\(r > 0.9\) for the pedigree examples provided).
Ancestral inbreeding coefficient
The ancestral inbreeding coefficient (\(F_{a}\), Ballou 1997) measures the probability of IBD of an individual for an allele that has been in homozygosity in at least one ancestor.
This parameter provides information not only about inbreeding, but can also be used to detect purging, since individuals with inbreeding \(F\) and ancestral inbreeding \(F_{a}\) are expected to be more fit than individuals with the same level of inbreeding but lower \(F_{a}\), given that the ancestors of the former have survived and reproduced despite their higher inbreeding (see Boakes & Wang 2005 and López-Cortegano et al. 2018 for analyses using this parameter).
Ancestral inbreeding can be estimated for an individual \(i\) with dam \(d\) and sire \(s\) as:
\[F_{a_{i}} = \frac{1}{2}[F_{a_{d}} + (1-F_{a_{d}})F_{d} + F_{a_{s}} + (1-F_{a_{s}})F_{s}]\] Alternatively, a gene-dropping simulation approach can be used, following Baumung et al. (2015), providing unbiased estimates of \(F_{a}\). This is because, above expression assumes that \(F\) and \(F_{a}\) are uncorrelated, which is not true.
Both approaches can be used with the function ip_Fa. Note that the computation of \(F_{a}\) requires estimating \(F\) in advance. Use argument Fcol to declare a column with \(F\) values if it has been computed and saved in advance (this will save time), or leave it blank to compute it on the go.
# F was pre-computed above
darwin %>%
purgeR::ip_Fa(Fcol = "Fi") %>%
dplyr::filter(names == "William Erasmus Darwin")
#> id dam sire names Fi Fa
#> 1 60 54 52 William Erasmus Darwin 0.06298828 0.001953125
# Compute F on the go (it won't be saved in the output)
# And enable genedropping
atlas %>%
purgeR::ip_Fa(genedrop = 1000, seed = 1234) %>%
dplyr::select(id, dam, sire, Fa) %>%
tail()
#> id dam sire Fa
#> 943 943 882 737 0.6475
#> 944 944 653 822 0.6000
#> 945 945 653 822 0.6000
#> 946 946 740 822 0.6050
#> 947 947 740 822 0.6075
#> 948 948 897 737 0.6340\(F_{a}\) can also be estimated based on population estimates of \(N_{e}\) and generation numbers, using the expression from López-Cortegano et al. (2018):
\[F_{a(t)} = 1 - (1-\frac{1}{2N})^{\frac{1}{2}t(t-1)}\] This can be achieved with the function exp_Fa (e.g. exp_Fa (Ne = 50, t = 50)).
Purged inbreeding coefficient
The purged inbreeding coefficient (\(g\)) gives the probability of IBD for deleterious recessive alleles. The reduction of \(g\) when compared to standard inbreeding depends on the magnitude of a purging coefficient (\(d\)) that measures the strength of the effective deleterious recessive component of the genome (García-Dorado 2012), so that \(d=0\) implies \(F=g\), and higher \(d\) (up to 0.5) means lower \(g\) in more inbred individuals. It can be calculated in pedigreed populations from the purged kinship coefficient (\(\gamma\)), in a similar way as standard inbreeding, following the methods described in García-Dorado (2012) and García-Dorado et al. (2016):
\[\gamma_{i,i} = \frac{1}{2}(1+g_{i})(1-2dF_{i})\]
\[\gamma_{i,j} = \frac{1}{2}(\gamma_{i,j_{d}}+\gamma_{i,j_{s}})(1-dF_{j})\] Where \(j_{d}\) and \(j_{s}\) are \(j\)’s mother and father respectively, and \(i\) is older than \(j\).
The function ip_g computes the purged inbreeding coefficient, given a value of \(d\). The choice of a proper value of \(d\) can however be complex. A separate vignette titled “Inbreeding and Purging Estimates” describes methods to help computing the inbreeding load as well as the purging coefficient.
atlas %>%
ip_F() %>%
ip_g(d = 0.48, Fcol = "Fi") %>%
dplyr::select(id, dam, sire, Fi, tidyselect::starts_with("g")) %>%
tail()
#> id dam sire Fi g0.48
#> 943 943 882 737 0.2350380 0.06066578
#> 944 944 653 822 0.2452226 0.08464775
#> 945 945 653 822 0.2452226 0.08464775
#> 946 946 740 822 0.2409467 0.07522799
#> 947 947 740 822 0.2409467 0.07522799
#> 948 948 897 737 0.2345642 0.06343757\(g\) can also be estimated based on population estimates of \(N_{e}\) and generation numbers, given a value of \(d\), using the expression from García-Dorado (2012):
\[g_{t} = [(1-\frac{1}{2N})g_{t-1}+\frac{1}{2N}](1-2dF_{t-1})\]
This can be achieved with the function exp_g (e.g. exp_g (Ne = 50, t = 50, d = 0.2)).
This is the last of functions related to inbreeding coefficients. Plotting together expected values of \(F\), \(F_{a}\) and \(g\) (assuming \(N_{e} = 25\) and an intermediate value \(d=0.25\)), the differences between the three coefficients become apparent.
data.frame(t = 0:50) %>%
dplyr::rowwise() %>%
dplyr::mutate(Fi = exp_F(Ne = 50, t),
Fa = exp_Fa(Ne = 50, t),
g = exp_g(Ne = 50, t, d = 0.25)) %>%
tidyr::pivot_longer(cols = c(Fi, Fa, g), names_to = "Type", values_to = "Inbreeding") %>%
ggplot(aes(x = t, y = Inbreeding, color = Type)) +
geom_line(size = 2) +
scale_x_continuous("Generations (t)") +
theme(legend.position = "bottom")
Opportunity of purging
Whereas previous purging methods focus on inbreeding measurements, opportunity of purging parameters calculate the potential reduction in individual inbreeding load, as a consequence of it having inbred ancestors (Gulisija and Crow 2007). The (total) opportunity of purging for an individual \(i\) (\(O_{i}\)) can be computed as:
\[O_{i} = \sum_{j}\sum_{k} (1/2)^{n-1} F_{j}\] Where \(j\) is every inbred ancestor of \(i\), \(k\) is every path from \(i\) to \(j\), and \(n\) is the number of individuals in the path (including \(i\) and \(j\)).
The expressed opportunity of purging depends on the probability of a given allele to be transmitted from an inbred ancestor \(j\) to \(i\), and thus on \(F_{i(j)}\). This is measured by the expressed opportunity of purging (\(O_{e}\)) as:
\[O_{ei} = \sum_{j} 2F_{i(j)} F_{j}\] In complex pedigrees (involving more than one inbred ancestor per path), these measures need to be corrected to discount the probability of purging measured in a close ancestor from that already calculated in a more distant ancestor. The function ip_op(complex=TRUE) does not perform Gulisija and Crow (2007) corrections, but instead applies an heuristic approach, only accounting for close ancestors \(j\), and ignoring contributions from far ancestors \(k\) such that \(F_{j(k) > 0}\).
The function ip_op can be used as:
arrui %>%
ip_op(Fcol = "Fi") %>%
dplyr::filter(target == 1) %>%
tidyr::pivot_longer(cols = c(Oe, Oe_raw)) %>%
ggplot() +
geom_point(aes(x = Fi, y = (value), fill = name), pch = 21, size = 3, alpha = 0.5) +
scale_y_continuous(expression(paste("Expressed opportunity of purging (", O[e], ")", sep=""))) +
scale_x_continuous("Inbreeding coefficient (F)") +
scale_fill_discrete("")
#> Computing partial kinship matrix. This may take a while.
The plot shows the increase of the expressed opportunity of purging with the inbreeding coefficient for individuals in the reference population (last two cohorts). For individuals with the lowest inbreeding values, the corrected (\(O_{e}\)) and uncorrected (\(O_{e \_raw}\)) have the same value, but as time progresses \(O_{e \_raw}\) becomes larger than \(O_{e}\). Both values are useful when determining potential reduction in the individual inbreeding load (see more exhaustive examples in López-Cortegano 2022, and the ‘Inbreeding and Purging Estimates’ vignette).
