true vs. estimated breeding values for bias estimation

[alanmaxsp]

Dear Chris and Gregor,

I am running a simulation that aims to mimic a beef cattle breeding program under additive + dominance effects. To preliminarily test the whole simulation I first simulated only 1 trait with additive effects (no dominance).

I have two questions related to true genetic values in AlphaSimR:

-In the preliminary analyzes I made an estimate of bias (both estimated and true) and I suspect that there may be some scaling (or means) problem between the EBVs and the TBVs (pop@gv); Is this possible? Should I do some kind of weighting or correction to the true genetic values?. The simple model I used in blupf90 is: y = overral mean + a + e.
I understand that to calculate accuracies (cor(EBV,TBV)) the scale is not a problem, but in the case of bias it is.

-Depending on the answer above, I'd like to get a thought from you about the A+D scenario.

Thanks a lot!
Alan

[gregorgorjanc]

@alanmaxsp one thing to note is that in AlphaSimR gv contain the base population trait mean (you are adding this when you are specifying the trait). Is that the “offset” that you might be alluding too?

Also note that GV are genetic values, not breeding values, but of course under additive these are the same in the first generation. Going onwards with generations AlphaSimR returns breeding values as they are defined (deviations from the population at hand) and as such cannot be compared between population objects - unless you merge the objects) - many animal breeders will be confused by this, hence I am highlighting it. Happy to say more.

[alanmaxsp]

Gregor, thank you very much for your comments.
Could you please check the simple example below?. Happy to hear your thoughts.
thanks in advance!

`

check how Alphasimr simulates genetic values

#Create founder haplotypes
founderPop = quickHaplo(nInd=10, nChr=1, segSites=10)

#Set simulation parameters
SP = SimParam$new(founderPop)
SP$addTraitAD(10, mean = 1, meanDD = 0, varDD = 0)

#Create a "fake" base population (note mean genetic values=1)
pop = newPop(founderPop, simParam=SP)

#checking means
SP$traits[[1]]@intercept; mean(genParam(pop)$gv)

#true genetic values:
pop@gv

#true additive effects in the base population (pop):
genParam(pop)$gv - SP$traits[[1]]@intercept
genParam(pop)$gv_a

#true breeding values in the base population (pop):
(genParam(pop)$gv_a + SP$traits[[1]]@intercept) - mean(genParam(pop)$gv_a + SP$traits[[1]]@intercept)
bv(pop)
gv(pop)-mean(gv(pop))

################################# PLEASE CHECK IT #####################
#Create 2nd population (or gen) derived from the pop:
pop2 = selectCross(pop, nInd=5, nCrosses=50, simParam=SP, use = "gv")

#true additive effects in pop2:
gv(pop2) - SP$traits[[1]]@intercept
genParam(pop2)$gv_a

#true genetic values in pop2:
genParam(pop2)$gv_a + SP$traits[[1]]@intercept
gv(pop2)

#true breeding value in pop2 relative to current population (pop2):
gv(pop2) - mean(gv(pop2))
bv(pop2)

#true breeding values in pop2 relative to both populations (pop y pop2):
pops <- c(pop,pop2)
gv(pops)[11:60] - mean(gv(pops))
bv(pops)[11:60]

#true breeding values in pop2 relative to fake base population (pop)???:
(genParam(pops)$gv_a[11:60] + SP$traits[[1]]@intercept) - mean(genParam(pops)$gv_a[1:10] + SP$traits[[1]]@intercept)
gv(pops)[11:60] - mean(gv(pops)[1:10])

#if the above is true
bv_mod <- function (pop, baseref = NULL,simParam = NULL)
{
if(!is.null(baseref)){
genParam(pop, simParam = simParam)$gv - mean(genParam(baseref, simParam = simParam)$gv)
}else{
genParam(pop, simParam = simParam)$bv
}
}
bv_mod(pop = pop2, baseref = pop)
bv_mod(pop = pop2, baseref = NULL)`

[gregorgorjanc]

@alanmaxsp I don't know what exactly I should be checking ;) You will have to be a bit more explicit and show the code and outputs to bring home what you are puzzled with - yes, it takes time, but it's worth it;)

[alanmaxsp]

Sorry for not being too explicit in the (very long) example!
Specifically, I would need to know if it is correct to get the true breeding values relative to some other base population in AlphaSimR like this:

#genetic values in a fake base population
gv(pop)
            [,1]
 [1,]  1.1418692
 [2,]  0.9108222
 [3,]  0.8471356
 [4,]  0.8065713
 [5,]  1.7588564
 [6,] -0.8098058
 [7,]  2.5910585
 [8,] -0.4806213
 [9,]  2.1530853
[10,]  1.0810286
mean_basepop = mean(gv(pop))
[1] 1

#genetic values in the population target (derived from the fake base population)
gv(pop2)
          [,1]
 [1,] 2.270941
 [2,] 1.126645
 [3,] 2.408194
 [4,] 1.081029
 [5,] 3.105415
 [6,] 2.895279
 [7,] 1.023082
 [8,] 1.824111
 [9,] 2.901191
[10,] 1.614032

#Get true breeding values in pop2 relative to fake base population (pop)
gv(pop2) - mean_basepop
            [,1]
 [1,] 1.27094094
 [2,] 0.12664452
 [3,] 1.40819419
 [4,] 0.08102860
 [5,] 2.10541474
 [6,] 1.89527871
 [7,] 0.02308172
 [8,] 0.82411058
 [9,] 1.90119093
[10,] 0.61403151

[gaynorr]

Hi Alan,

What you need here is the average effect of an allele substitution (alpha), so that you can reconstruct the breeding values manually. The most common formula provided in textbooks makes the assumption of Hardy-Weinberg equilibrium (HWE): alpha = a + d * (q - p). Where a is the additive effect, d is the dominance effect, p is the frequency of the 1 allele, and q is the frequency of the 0 allele. Without assuming HWE, the formula is alpha = a + d * (q - p) * (1 - F) / (1 + F). Where F is the coefficient of inbreeding. You'll need this second formula to reconstruct the value used by AlphaSimR.

You can then calculate the breeding values for the other population by extracting its QTL genotypes and subtracting 2p from their value before multiplying by alpha. The subtraction of 2p is critical here, because you're interested in bias. Here's an example below.

library(AlphaSimR)

# Create founder haplotypes
founderPop = quickHaplo(nInd=100, nChr=10, segSites=1000)

# Create a trait with additive and dominance effects
SP = SimParam$
  new(founderPop)$
  addTraitAD(1000, meanDD=0.2, varDD=0.1)

# Create an initial population from founder haplotypes
pop1 = newPop(founderPop)

# Create a second population by mating the first
pop2 = randCross(pop1, nCrosses=100)


# Extract QTL genotypes from pop1 to calculate p and q
qtlGeno = pullQtlGeno(pop1) 
p = colMeans(qtlGeno) / 2
q = 1 - p

# Calculate F using expected heterozygosity for HWE (2 * p * q)
expHet = 2*p*q # Expected heterozygote frequency.
obsHet = colMeans( 1 - abs(qtlGeno-1) ) # Observed heterozygote freq.
F = 1 - obsHet / expHet

# Calculate average effect of an allele subsitution
a = SP$traits[[1]]@addEff
d = SP$traits[[1]]@domEff
alpha = a + d * (q-p) * (1-F) / (1+F)

# Recreate breeding values for pop1
qtlGenoCenter = sweep(qtlGeno, 2, 2*p)
BV_pop1 = qtlGenoCenter %*% alpha

# Check validity of values
mean(BV_pop1) # Should be zero within numeric precision
cor(BV_pop1, bv(pop1)) # Should be 1

# Calculate values for pop2
qtlGeno2 = pullQtlGeno(pop2)
qtlGenoCenter2 = sweep(qtlGeno2, 2, 2*p)
BV_pop2 = qtlGenoCenter2 %*% alpha

[alanmaxsp]

Thanks for the excellent reply Chris!

So, with a trait without dominance effects, it is expected that TBV (calculated that way) == TGV (as reported by alphasimr), as long as the desired mean of the base population is zero; am I right?

Regards!

[gaynorr]

Yes, as long as your base population is the same population provided to SimParam to define the trait. When the trait is defined, an intercept is set to achieve your requested mean in this population. You can access this intercept directly for the first trait using this code: SP$traits[[1]]@intercept.