26
$\begingroup$

I've set up some dummy data in R which makes 40 genetically related lines, they are all siblings within a line so are genetically related by a factor of ½ thus additive genetic variance should be twice the variance explained by line. For the lines there are 200 individuals being measured, for three characters/traits. The first trait has low phenotypic variance, the second has high environmental variance, and the third has high genetic variance.

rm(list=ls())
re = 200 # replicate individuals per line
li = 40  # lines

# setup
set.seed(5)
data = data.frame(rep(1:li, each = re))
colnames(data)="Line"
library(nlme)
library(lme4)
par(mfrow=c(1,3))

# trait 1: little variance (Va or Ve)
data$Trait = rnorm(li*re,10,1)
boxplot(data$Trait~data$Line, ylim=c(0,20), main = expression("Low V"[A]*"& Low V"[E]))
var1 = var(data$Trait); mean1 = mean(data$Trait); m1 = lmer(data$Trait ~ (1|data$Line))

# trait 2: high enivronmental variance, little Va
data$Trait = rnorm(li*re,10,1); data$Trait = data$Trait + rnorm(re*li,0,3)
boxplot(data$Trait~data$Line, ylim=c(0,20),main = expression("Low V"[A]*"& High V"[E]))
var2 = var(data$Trait); mean2 = mean(data$Trait); m2 = lmer(data$Trait ~ (1|data$Line))

# trait 3: high additive genetic variance, little Ve
data$Trait = rnorm(li*re,10,1); data$Trait = data$Trait+rep(rnorm(li,0,3),each=re)
boxplot(data$Trait~data$Line, ylim=c(0,20), main = expression("High V"[A]*"& Low V"[E]))
var3 = var(data$Trait); mean3 = mean(data$Trait); m3 = lmer(data$Trait ~ (1|data$Line))

Plots of the three traits,

enter image description here

Then to estimate additive variance ($V_A$) I extract the line variance ($V_L$) from the lmer model and double it,

# line variances (variance in additive effect of each haploid genome)
m1_line = unlist(VarCorr(m1))[[1]]; 
m2_line = unlist(VarCorr(m2))[[1]]; 
m3_line = unlist(VarCorr(m3))[[1]]

# additive variance (double the line variance because it is a hemiclone)
m1_add = 2*m1_line; 
m2_add = 2*m2_line; 
m3_add = 2*m3_line

Residual variances should be (assuming perfect experimental design, no measurement error etc.) the estimate of environmental variance ($V_E$) and phenotypic variance ($V_P$) should be the sum of $V_L$ and $V_E$,

# residual variance
m1_res = attr(VarCorr(m1), "sc")^2
m2_res = attr(VarCorr(m2), "sc")^2
m3_res = attr(VarCorr(m3), "sc")^2

# phenotypic variance
m1_phe = m1_line + m1_res
m2_phe = m2_line + m2_res
m3_phe = m3_line + m3_res

Heritability is additive variance divided by the phenotypic variance

$h^2 = V_A / V_P$

But I think it is correct in this case to use line variance rather than additive variance (if someone could explain in an answer that would be useful), so I've done,

# heritability (line/ (line+ residual))
m1_h2 = m1_line/ m1_phe
m2_h2 = m2_line/ m2_phe
m3_h2 = m3_line/ m3_phe
m1_h2; m2_h2; m3_h2

My question(s):

Is it appropriate to use the lmer function in R to extract variance components in this manner?

Have I calculated $V_A$, $V_E$, $V_P$, $h^2$ correctly? I think $V_A$ and $V_E$ are correct, $V_P$ could perhaps be the sum of $V_A$ and $V_E$ rather than $V_L$, and subsequently $h^2$ may be $V_A/V_P$.

$\endgroup$
  • $\begingroup$ Can you add a flowchart or just methodology instead of/in addition to the code. $\endgroup$ – WYSIWYG Oct 9 '14 at 9:47
  • 3
    $\begingroup$ I deleted my answer. I am too unsure of what I was saying! Sorry. I might give another try thinking about that later when I'll have the opportunity to use R. $\endgroup$ – Remi.b Oct 9 '14 at 15:05
0
$\begingroup$

From what I understood of your code and what you are asking I am guessing that you do the following:

Generating a virtual set of 40 individuals (lines) of which you have 200 measurements (repetitions). You say that they are full siblings, so they share both parents. Then you use the lmer function (which I am not familiar with) to give you the total variance, within-group variance and between-group variance, ($V_{T}, V_W, V_B$ respectively). What you call $V_L$ would be the $V_W$.

We know that $V_P=V_A+V_D+V_E$

In full siblings we also know that $V(A,A')=\frac{1}2V_A+\frac{1}4V_D$ where $A$ is the additive value and $D$ is the dominance value because they share a quarter of their gene combinations.

If you are treating with full siblings the $V_W$ should be then the $V_E$ and the $V_B$ would be the $\frac{1}2V_A+\frac{1}4V_D$ as you are seeing the variance between siblings.

So no, you are not calculating correctly the $V_A$ because you are missing your $V_D$ on your counting. The code itself looks fine for what I understood, but maybe you should take that to corssvalidated with the R tag

$\endgroup$
  • $\begingroup$ Thanks - I actually use Hemiclone lines, which are like full sibs (related by a half) but have no dominance effect, I just wrote siblings rather than having to explain hemiclones. If then half sibs were used, the code I use would correct except for $V_A = 4V_L$ right? And likewise if there was no dominance then for full sibs $V_A = 2V_L$ $\endgroup$ – rg255 Aug 21 '15 at 8:08
  • $\begingroup$ I could check it thoroughly but my gut says that $V_A=2V_L$ in this case as there would be no dominance effect then, but in the rest of the effects it would be full sibs. $\endgroup$ – Athe Aug 21 '15 at 8:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.