부모-자식 데이터를 통한 유전율 추정

부모-자식 데이터를 통한 유전율 추정

본 포스팅은 부모-자식 데이터를 통해 heritability를 추정하는 것을 포스팅한다. 아래 데이터는 parent와 child의 height를 나타내는 수치이다. 옥수수의 키라고 생각하자. 알고 싶은 것은 child의 키가 유전적 요소로 설명되는 정도이다. 우리는 이 수치를 heritability(h^2)라는 값으로 표현한다. 데이터는 이곳에서 다운로드 할 수 있다.

\[ h^2 = \frac{V(G)}{V(Y)}\]

library(data.table)

## Warning: package 'data.table' was built under R version 3.4.3

load("C:\\workspace\\r\\genetics\\data\\mod7_data.R")
  head(Galton)

##     parent child
  ## 1 71.40788  61.7
  ## 2 68.57945  61.7
  ## 3 64.33681  61.7
  ## 4 62.92260  61.7
  ## 5 62.21549  61.7
  ## 6 67.16524  62.2

데이터에 약간의 random error를 주어 아래의 그래프를 통해 표현해보자.

plot(Galton$parent+rnorm(928,0,.1),Galton$child+rnorm(928,0,.1),xlab="Father's Height",ylab="Child's Height")
  abline(v=mean(Galton[,1]),col=2,lty=3,lwd=3)
  abline(h=mean(Galton[,2]),col=2,lty=3,lwd=3)
  abline(a=0,b=1,col=2,lty=3,lwd=3)

summary(lm(Galton$child~Galton$parent))

## 
  ## Call:
  ## lm(formula = Galton$child ~ Galton$parent)
  ## 
  ## Residuals:
  ##     Min      1Q  Median      3Q     Max 
  ## -7.8050 -1.3661  0.0487  1.6339  5.9264 
  ## 
  ## Coefficients:
  ##               Estimate Std. Error t value Pr(>|t|)    
  ## (Intercept)   36.87187    1.98827   18.55   <2e-16 ***
  ## Galton$parent  0.45700    0.02909   15.71   <2e-16 ***
  ## ---
  ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  ## 
  ## Residual standard error: 2.239 on 926 degrees of freedom
  ## Multiple R-squared:  0.2105, Adjusted R-squared:  0.2096 
  ## F-statistic: 246.8 on 1 and 926 DF,  p-value: < 2.2e-16

회귀계수가 0.457이다. 이 때, 유전율은 0.457*2 = 약 92%로 추정된다. 왜 그런지 살펴보기 위해 우선 beta의 추정치를 아래와 같이 회귀식의 양변에 Cov(, X)를 취해 구해보자.

\[ Y = \alpha + \beta X \] \[ Cov(Y, X) = \beta Cov(X,X) \] \[ \hat{\beta} = \frac{Cov(Y,X)}{Var(X)} \]

Phenotype Y에 대해 아래와 같은 일반적인 Phenotype에 대한 모델을 세우자. 이 때, A,D,C,E 각각의 컴포넌트들에 대한 것은 나중에 살펴보고 i, j 두 개체의 Covariance만 살펴본다. A : additive effect D : dominance effect C : shared environment effect E : residuals \[ Y = A + D + C + E\] r은 relatedness, k2는 genome상의 IBD2의 비율, a는 shared environment를 나타내는 수치이다. 만약 shared environment이면 a=1, 아니면 a=0이다. 

이 Covariance 모델은 어디까지나 가정이다. 이 가정 하에서 각각의 Variance를 추정하고 이를 통해 heritability를 추정하게 된다.

\[ E[Cov(Y_i, Y_j)] = rVar(A) + k_2Var(D) + aVar(C)\]

이제 현재 데이터에서 위 식을 적용하자. 이 때, r=0.5, k2=0, a=0을 가정하자. (부모-자식 간의 IBD2=0이므로 k2=0이다. ) 그러면 랜덤으로 i번째 부모-자식 pair를 골랐을 때, Covariance 모델은 아래와 같다. (P : Parent, C : Child의 약자)

이러면 Var(D)가 사라지기 때문에 우리는 쉽게 Var(A)를 데이터를 통해 추정할 수 있다.

\[ E[Cov(P_i, C_i)] = 0.5Var(A)\] \[ Var(A) = 2E[Cov(P_i, C_i)] \]

print(var(Galton$parent))

## [1] 6.389121

print(var(Galton$parent, Galton$child)) # Sample Covariance

## [1] 2.919806

따라서 위 데이터에는 D(Dominance effect가 없으므로 아래와 같이 유전율을 추정할 수 있다.) \[ h^2 = \frac{V(G)}{V(Y)} = \frac{V(A)}{V(P)} \]

따라서 2.91*2 = 5.82 = Var(A) 이다. Var(P) = 6.39 (Phenotype의 분산을 V(P)로 추정) 이므로 이 모델하에서 추정한 유전율은 아래와 같다.

\[ \frac{V(A)}{V(P)} = 5.82/6.39 = 0.91 \]

print(paste0("Heritability : ", 2*var(Galton$parent, Galton$child) / var(Galton$parent)))

## [1] "Heritability : 0.913992906289406"

따라서 부모-자식 데이터에 이러한 모델을 통해 유전율을 추정할 때 유전율은 회귀계수의 2배를 해주면 된다.

참고

https://jvanderw.une.edu.au/AGSCcourse.htm