OPEN ACCESS JOURNAL
Peer Reviewed
eISSN
선형모형을 이용하여 단일염기변이(Single Nucleotide Polymorphism, SNP)의 효과를 추정하는 방법론은 유전체 정보를 이용하여 개체의 유전능력을 추정하는 가축육종(Animal Breeding) 분야에서 매우 중요한 내용일 것이다. 이번 장에서는 선형모형의 일종인 회귀분석과 분산분석의 개념적 차이 및 개별 방법론에 대한 중요한 점에 대해서 학습한다. 회귀분석과 분산분석(ANOVA)은 모두 종속변수와 독립변수 사이의 관계를 분석하는 통계기법이지만, 접근 방식에 약간의 차이가 있다. 회귀분석은 보통 독립변수가 연속형(continuous variable)일 때 선형적 관계(회귀식의 기울기)를 추정하는 데 초점을 두고, 분산분석은 독립변수가 범주형(categorical variable)일 때 집단 간 평균의 차이에 초점을 둔다. 예를 들어, SNP 유전자형이 연속형 코딩(예: 0,1,2)으로 주어지면 회귀분석을 통해 종속변수에 대한 기울기(allele 당 효과)를 추정할 수 있고, SNP 유전자형이 범주형 요인(AA, Aa, aa의 3수준)으로 주어지면 분산분석을 통해 세 집단 평균 차이에 대한 F-검정을 수행하게 된다.
두 방법은 모델의 기본 구조가 동일한 일반선형모델(GLM)이며, 독립변수를 어떻게 취급하는가에 따라 출력 해석이 달라진다. 회귀분석에서는 회귀계수(예: SNP 대립유전자 하나 증가당 효과)를 해석하고, 분산분석에서는 집단별 평균과 집단 간 분산 대비 집단 내 분산의 비율(F값)을 해석한다. 결과적으로 회귀분석과 일원분산분석은 동일한 모형을 다른 형태로 표현한 것이며, 회귀분석은 개별 계수(예: 더미변수나 기울기)에 대한 t-검정을 제공하고, 분산분석은 전체 요인의 효과에 대한 F-검정을 제공한다. 요약하면, “ANOVA는 모든 예측변수가 범주형인 회귀분석의 특수한 경우”라고 볼 수 있다.
SNP 마커 분석 적용: SNP의 유전자형을 어떻게 모델링하는가에 따라 두 접근이 모두 활용된다. 예를 들어 SNP 유전자형을 0/1/2로 (대립유전자 개수에 따라) 코딩하면 회귀분석으로 상가적 효과(additive genetic effect)를 추정할 수 있고, SNP 유전자형을 AA/Aa/aa 세 그룹으로 취급하면 분산분석으로 그룹 간 평균의 유의성을 평가할 수 있다. 두 방법은 동일한 데이터를 다루며 결과도 밀접히 관련되어 있다. 실제로 회귀분석으로 구한 더미변수(dummy variable) 회귀계수로부터 그룹 간 평균차를 구할 수 있고, 분산분석의 F-검정은 해당 회귀모형에서 더미변수(dummy variable) 계수들의 공동 유의성 검정과 동일하다. 차이는 연구자가 관심을 가지는 표현 방식에 있는데, 집단 평균 차이에 직접 관심이 있다면 ANOVA의 F-검정과 사후검정을 주로 보고, 대립유전자 하나당 효과 크기에 관심이 있다면 회귀계수(기울기)를 보고 해석하는 식이다.
SNP 마커의 유전자형이 AA, Aa, aa로 주어졌을 때, 이를 세 그룹으로 보고 일원분산분석(one-way ANOVA)을 수행하는 절차는 다음과 같다.
* 데이터 준비: 각 개체의 유전자형을 범주형 변수(factor)로 지정한다. 예를 들어, 유전자형이 AA, Aa, aa 중 하나인 열을 인코딩한다. 필요시 Aa와 aA처럼 표현이 섞여 있다면 동일하게 취급하도록 정리한다. (두 표현은 생물학적으로 동일한 이형접합).
* 집단별 기술통계: 분산분석 전에 각 유전자형 그룹의 표본크기(N)와 평균, 분산 등을 산출한다. 이를 통해 그룹 간 평균 차이가 어느 정도인지, 분산이 유사한지 등을 파악할 수 있다. (ANOVA의 기본 가정 중 하나는 각 그룹의 분산이 동일하다는 등분산 가정)
* 분산분석 실행: 종속변수 ~ 유전자형을 모형으로 하여 일원분산분석을 수행한다. 예를 들어 R에서는 aov(phenotype ~ genotype) 또는 lm(phenotype ~ genotype)를 사용하는데, 이때 genotype 변수는 세 수준(AA, Aa, aa)의 팩터(factor)로 취급한다.
* 분산분석 F-검정: ANOVA 결과에서는 집단간 제곱합과 오차 제곱합을 비교하는 분산비(F)가 계산된다. 분산분석은 기본적으로 집단간 변동이 무작위 오차변동에 비해 충분히 큰지를 검정하며, 유전자형에 따른 집단 평균 차이가 통계적으로 유의미한지 p값을 확인한다. 예를 들어, 아래 ANOVA 결과에서 유전자형 효과에 대한 F-통계량과 p값을 해석한다.
| 분산 | 제곱합 | 자유도 | 평균제곱 | F검정통계량 | 유의확률 |
|---|---|---|---|---|---|
| 집단간 | 17.69 | 2 | 8.85 | 119.01 | 0.000 |
| 집단내 | 7.21 | 97 | 0.07 | ||
| Total | 24.90 | 99 |
위 예시에서 유전자형 요인은 2개의 자유도(df=2)를 가지며, F값은 약 119.01, p값은 0.000으로 나와 있다. 이는 유전자형 간 평균 차이가 유의미함을 보인다. 분산분석의 귀무가설(H0)은 “세 가지 유전자형 집단의 모평균은 모두 같다”이며, 대립가설은 “적어도 하나의 집단 평균이 다르다”로 표현한다. 따라서 p값이 충분히 작으면 (예: < 0.05) 유전자형에 따라 그룹 평균에 차이가 있다고 결론 내릴 수 있다.
* 사후검정 및 효과추정: ANOVA의 F-검정이 유의한 경우, 어떤 그룹 간 차이가 유의한지를 추가로 확인하기 위해 Tukey HSD 등의 사후검정을 수행한다. 또한 집단별 평균값의 차이 자체가 마커의 효과 추정치가 되므로, AA 대비 aa의 평균 차이 (또는 기준 유전자형 대비 다른 유전자형의 차이)를 구해 마커 효과의 크기를 볼 수 있다.
요약하면, 유전자형을 세 그룹으로 나누어 ANOVA를 수행하는 것은 각 유전자형별 표현형 평균을 비교하는 것이며, 이는 해당 SNP가 표현형 변이에 영향을 주는지 탐색하는 기본적인 방법이다. 이때 ANOVA에서는 개체 간의 총 변동을 유전자형 집단 간 변동과 집단 내(오차) 변동으로 분해하여, 집단간 변동이 얼마나 큰 비율을 차지하는지를 F값으로 평가한다. 유전자형에 따른 효과가 크다면 집단간 평균차가 커져 F값이 커지고, 효과가 없다면 세 평균이 유사하여 F값이 작게 나온다.
상가적 유전모형(additive genetic model)이란 대립유전자의 효과가 서로 대등하게 표현형에 기여하여, 이형접합의 표현형이 두 동형접합의 중간 정도로 나타나는 모형을 말한다. 즉 한 대립유전자씩의 기여분이 균등하다고 가정하는 모델이다. 이러한 모형에서는 마커의 효과를 대립유전자 한 개가 표현형에 주는 증가분으로 정의할 수 있다. 실질적으로는 특정 대립유전자(주로 minor alleles, 희귀 대립유전자)를 기준으로 0, 1, 2개 보유한 경우를 연속등급으로 놓고 회귀계수로 효과를 추정한다. 예를 들어 어떤 SNP의 두 대립유전자를 A(주류 allele)와 a(희귀 allele)라고 할 때, a allele을 하나도 가지지 않은 유전자형을 0, 이형접합(Aa)을 1, aa 동형접합을 2로 코딩하여 회귀분석을 하면 기울기(회귀계수)가 대립유전자(allele) a 한 개당 평균적인 표현형 변화량을 추정하게 된다. 이를 통상 대립유전자 치환 효과(allelic substitution effect) 또는 상가적 효과라고 부른다.
이론적으로 상가적 효과는 두 동형접합의 평균 차이의 절반으로 정의되는데, 표현형 값의 모집단 평균이 동형 유전자형$(M)_{AA}$, 동형 유전자형$(m)_{aa}
$ (각각 AA, aa 유전자형의 평균)이고 이형접합 유전자형_Aa의 평균이 $m_{Aa}$라고 하면, 상가적 효과를 $a=\dfrac{m_{aa}-m_{AA}}{2}$로 정의할 수 있다. 이는 allele a를 하나 더 가질 때 기대되는 증가분을 의미한다. 완전 상가적 모형이라면 이형접합의 평균은 정확히 두 동형접합 평균의 중간값이 되며 , $m_{Aa}\approx\dfrac{m_{AA}-m_{aa}}{2}$ 오차 범위 내에서 $m_{Aa}\approx m_{AA}+a\approx m_{aa}-a$가 성립한다. 이 경우 a allele의 효과는 A allele의 효과와 크기는 같고 방향은 반대이며, 표현형에 대한 aa 유전자형의 총 효과는 2a (AA 대비)로 볼 수 있다.
만약 우성(dominance) 효과가 존재한다면, 이형접합의 표현형이 완전히 중간이 아니라 한쪽 동형접합에 더 가깝거나 멀어지게 된다. 이를 수량적으로 표현하면 우성효과의 편차(dominance deviation) d를 $m_{Aa}$와 중간값의 차이로 정의할 수 있다. $d=m_{Aa}-\dfrac{m_{AA}+m_{aa}}{2}$. 순수 상가적 모형에서는 $d \approx 0$이며, $m_{Aa}$가 정확히 중간값이 된다. 우성의 효과가 있다면 $d \ne 0$이고, 양수이면 이형접합이 aa쪽에 가깝고 음수이면 AA쪽에 가깝다는 뜻이다. 상가적 모형 가정 하에서는 d=0으로 보고, 마커 효과는 전적으로 상가적 효과 a로 설명된다.
이러한 개념을 통해 마커 효과의 추정이 이루어진다. 상가적 효과 $a$는 위 정의대로 $AA$와 $aa$의 평균차의 절반으로 계산할 수 있고, 이는 회귀분석에서 구한 대립유전자 1개당 효과와 동일하다. 예를 들어 앞의 SNP에서 allele a를 0,1,2로 코딩하여 회귀분석을 했다면 회귀계수 $\hat{\beta}$가 바로 추정된 상가적 효과이며, 이는 (표본에서의 $\dfrac{\overline{m_{aa}}-\overline{m_{AA}}}{2}$와 거의 같게 계산된다. 상가적 효과는 흔히 육종에서의 유전력 계산이나 BLUP 산출 등에 중요하며, 유전형의 Breeding value로 연결되는 개념이다.
정리하면, 상가적 유전모델에서는 대립유전자 한 개의 효과를 일정한 값으로 간주하며, SNP 마커 효과는 “해당 SNP의 대립유전자가 한 개 늘어날 때 기대되는 표현형 변화”로 정의된다. 이를 위해 SNP 유전자형을 0/1/2로 코딩하여 회귀분석을 수행하며, 이렇게 추정된 회귀계수는 곧 마커의 상가적 효과 추정치가 된다. 만약 이 값이 통계적으로 유의하면 해당 SNP의 효과가 유의미하게 존재한다고 해석한다. 반대로 분산분석에서는 한 번에 세 유전자형 평균을 비교하므로, 상가적 + 우성의 효과를 합쳐 유전자형 전체의 영향을 검정하게 된다. 상가적 모형이 적합하다면 회귀분석(1 df)이 분산분석(2 df)보다 검정력이 높을 수 있고, 우성의 효과를 무시할 수 없으면 분산분석이 더 포괄적인 검정을 제공한다.
유전자형을 범주형과 연속형으로 간주한 간단한 예제를 이용하여 ANOVA모델과 Regression 모델의 일차방정식의 행렬식을 구현하고 모델이 어떻게 마커효과의 해(solution)을 추정하는지 학습한다. 각 방법의 이론적 배경, 행렬식 구현, 실습 예제 이용하여 설명한다.
일반 선형모형(GLM)은 여러 요인($X$)이 관찰값($Y$)에 미치는 영향을 설명하기 위한 대표적인 통계모형이다. 행렬식을 통해 일반 선형모형을 표현하면 아래와 같다.
$$\mathbf{Y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon}$$
$\mathbf{Y} (n \times 1)$: 관찰된 표현형 벡터 (예: 개체별 체중)
$\mathbf{X} (n \times p)$: 디자인 행렬(어떤 개체가 어떤 요인에 속하는지 정리한 표)
$\boldsymbol{\beta} (p \times 1)$: 고정효과 벡터(요인별 효과 크기)
$\boldsymbol{\varepsilon}(n \times 1)$: 오차 벡터, $E[\boldsymbol{\varepsilon}]=0$, $Var[\boldsymbol{\varepsilon}]=\sigma^2 I_n$
이 식을 풀어 $\boldsymbol{\beta}$를 추정하려면, ‘다음의 추정방정식(estimation equations)’을 사용한다. 여기서는 아래의 $\boldsymbol{\beta}$추정 방정식을 유도하는 것은 생략하도록 한다.
$$\mathbf{X}^T \mathbf{X} \hat{\boldsymbol{\beta}}=\mathbf{X}^T \mathbf{Y} \Rightarrow \hat{\boldsymbol{\beta}}=(\mathbf{X}^T \mathbf{X})^{-1}\mathbf{X}^T \mathbf{Y}$$
위의 행렬식을 이용하여 $\mathbf{X}$행렬과, $\mathbf{Y}$ 행렬식을 이용하여 $\boldsymbol{\beta}$를 추정하는 행렬식 계산 과정(곱셈, 전치, 역행렬 등)은 뒤 예제에서 구체적으로 살펴보도록 한다.
| 개체 | 체중(kg) | 유전자형 | 대립유전자 수(G) |
|---|---|---|---|
| 1 | 22 | AA | 0 |
| 2 | 24 | AA | 0 |
| 3 | 26 | AA | 0 |
| 4 | 28 | AG | 1 |
| 5 | 30 | AG | 1 |
| 6 | 32 | AG | 1 |
| 7 | 39 | GG | 2 |
| 8 | 38 | GG | 2 |
| 9 | 40 | GG | 2 |
본 예제의 대상 개체 수는 총 9마리($n=9$)이다. 처음 3마리는 AA유전자형, 그리고 그 다음 3마리는 AG유전자형 나머지 3마리는 GG유전자형을 가지고 있다. 평균 체중은 $\bar{Y}=31kg$이다. 이 간단한 데이터를 이용하여 두 모델이 어떻게 $X$ 행렬을 구성하고 $\beta$를 추정하는지 살펴봅니다.
이론적 배경
ANOVA(분산분석)는 집단별 평균 비교에 최적화된 방법이다. 이론적으로는 ‘AA’, ‘AG’, ‘GG’ 세 개 그룹의 평균 간 차이를 검정하는 것으로 볼 수 있다. 모형식은 아래와 같다.
모형식:
$$y_i=\beta_0 + \beta_1 I(AG)_i + \beta_2 I(GG)_i+e_i$$
$I(AG)_i$, $I(GG)_i$: AG, GG이면 1, 아니면 0이 되는 지시 변수(dummy)
$\beta_0$: AA 그룹 평균
$\beta_1$: AG vs AA 평균 차이
$\beta_2$: GG vs AA 평균 차이
이 예제를 이용해서 이제 행렬식을 구현해 보자. 유전자형을 설명하는 디자인 행렬 $\mathbf{X}_{ANOVA}(9 \times 3)$: 총 개체수 9마리와 유전자형 3개 따라서, 9개의 행과 3개의 열로 구성되어 있다. 그리고 우리가 구하고자 하는 해는 총 3개로 ANOVA 로 표시할 수 있다.
$$\mathbf{X}_{ANOVA}=\begin{pmatrix}
1 & 0 & 0 \\
1 & 0 & 0 \\
1 & 0 & 0 \\
1 & 1 & 0 \\
1 & 1 & 0 \\
1 & 1 & 0 \\
1 & 0 & 1 \\
1 & 0 & 1 \\
1 & 0 & 1
\end{pmatrix}, \boldsymbol{\beta}_{ANOVA}=\begin{pmatrix}
\beta_{0} \\ \beta_{1} \\ \beta_{2}
\end{pmatrix}$$
위의 $\mathbf{X}$ 행렬식과, 표현형인 $\mathbf{Y}$ 행렬식을 이용하여, 아래와 같은 전치행렬과 행렬의 곱을 이용하고, $\mathbf{X}$전치행렬과 $\mathbf{X}$행렬을 곱한 $\mathbf{X}^{T}\mathbf{X}$의 역행렬을 구하여 최종 $\boldsymbol{\beta}_{ANOVA}$를 구한다.
$$\mathbf{X}^{T}\mathbf{X}=\begin{pmatrix}
9 & 3 & 3 \\
3 & 3 & 0 \\
3 & 0 & 3
\end{pmatrix}, \mathbf{X}^{T}\mathbf{Y}=\begin{pmatrix}
279 \\ 90 \\ 117
\end{pmatrix}$$
이를 역행렬로 풀면:
$$\hat{\boldsymbol{\beta}}_{ANOVA}=(\mathbf{X}^{T}\mathbf{X})^{-1}\mathbf{X}^{T}\mathbf{Y}=\begin{pmatrix}
24 \\ 6 \\ 15
\end{pmatrix}$$
결과적으로, AA=24kg, AG=24+6=30kg, GG=24+15=39kg로 계산되고, 이를 직관적으로 설명하면, AA 그룹 세 개 체중 평균을 직접 구하면 24kg, AG 그룹은 AA 평균(24)에 +6을 더해, 30kg, GG 그룹은 AA 평균(24)에 +15를 더해, 39kg이다. 이를 활용하여 상가적 유전효과를 계산하면, (39kg-24kg)/2로 약 7.5kg으로 계산된다. 그러면 이제부터 회귀모형을 이용하여 해를 구해보자.
이론적 배경
Regression(회귀분석)은 상가적(additive) 효과를 가정하여, ‘대립유전자 개수’가 연속적으로 체중에 영향을 준다고 봅니다. 그리고 모형식은 아래와 같습니다.
모형식:
$$y_i = \beta_0+\beta_1 G_i+e_i$$
$G_i \in \{0,1,2\}$: AA=0, AG=1, GG=2
$\beta_0$: AA 평균
$\beta_1$: a 대립유전자 하나당 체중 변화량
이 자료를 이용하여 행렬식을 구현하면 다음과 같습니다. 먼저 디자인 행렬은 9개의 행과 2개의 열로 구성된다.
디자인 행렬 $\mathbf{X}_{REG} (9 \times 2)$ , 그리고 추정해야 하는 해는 2개이다.
$$\mathbf{X}_{REG}=\begin{pmatrix}
1 & 0 \\
1 & 0 \\
1 & 0 \\
1 & 1 \\
1 & 1 \\
1 & 1 \\
1 & 2 \\
1 & 2 \\
1 & 2
\end{pmatrix}, \boldsymbol{\beta}_{REG}=\begin{pmatrix}
\beta_0 \\
\beta_1
\end{pmatrix}$$
이 행렬을 이용하여 $\mathbf{X}$전치행렬, $\mathbf{X}$의 곱과, $\mathbf{X}$의 전치행렬과 $\mathbf{Y}$를 곱하여 아래와 같은 행렬을 구성하고,
$$\mathbf{X}^T \mathbf{X} = \begin{pmatrix}
9 & 9\\
9 & 15
\end{pmatrix}, \mathbf{X}^T \mathbf{Y} = \begin{pmatrix}
279\\
324
\end{pmatrix}$$
역행렬을 이용하여 해를 구하면 다음과 같다.
$$\hat{\boldsymbol{\beta}}_{REG} = (\mathbf{X}^T \mathbf{X})^{-1}\mathbf{X}^T \mathbf{Y} = \dfrac{1}{54} \begin{pmatrix}
15 & -9\\
-9 & 9
\end{pmatrix} \begin{pmatrix}
279\\
324
\end{pmatrix} = \begin{pmatrix}
23.5\\
7.5
\end{pmatrix}$$
결과적으로 AA≈23.5kg, a 한 개당 +7.5kg → AG(1)≈31kg, GG(2)≈38.5kg로 해가 ANOVA모델과 동일하게 계산되는 것을 확인할 수 있다. 즉, 절편 23.5는 AA그룹의 평균값, 그리고 기울기 7.5는 a 한 개를 더 가질 때마다 체중이 7.5kg 증가함을 의미한다.
| 구분 | ANOVA (범주형) | Regression (연속형) |
|---|---|---|
| 모델 가정 | 집단별 평균 차이 검정 | 대립유전자 수당 선형 증가 가정 |
| 모수(parameter) 수 | 3 ($\beta_{0}$, $\beta_{1}$, $\beta_{2}$) | 2 ($\beta_{0}$, $\beta_{1}$) |
| 해석 | 그룹별 절대평균 | 계수당 상대효과 |
| 검정력 | 우성효과(dominance) 포착 가능 | 상가적(additive) 효과 강하게 검출 |
| 확장성 | 범주형 요인 추가 용이 | 다변량 연속형 예측변수 확장 용이 |
ANOVA는 ‘어느 그룹이 얼마나 다른가’를 명확히 비교할 때, 우성효과(Dominance)나 비선형 효과를 탐색할 때 적합하다. Regression은 ‘유전자가 하나씩 더해질 때’라는 명확한 가정을 바탕으로, 여러 마커를 동시에 분석하거나 예측 모델에 사용할 때 유용하다. 두 방법 모두 $\mathbf{X}^{T} \mathbf{X} \cdot \beta = \mathbf{X}^{T}\mathbf{Y}$ 행렬 계산으로 통일되므로, 디자인 행렬 구성만 달리하면 쉽게 전환 가능하다. 이 예제를 통해, 동일한 데이터라도 디자인 행렬 $\mathbf{X}$의 구조를 바꿈으로써 서로 다른 관점(집단 평균 vs 연속적 기울기)으로 분석할 수 있음을 보여주었다.
이 강의자료에서는 약 100마리 한우의 근내지방함량(D_IMF) 데이터를 활용하여 SNP 마커(유전자형)의 효과를 추정하는 두 가지 주요 통계적 방법ANOVA(분산분석)과 Regression(회귀분석)을 비교한다. 실제 한우 등심의 근내지방함량(D_IMF)를 표현형(Y)으로 설정하고, 유전자형 그룹 간 평균 차이를 검정하고, 상가적(additive) 모형 하에서 대립유전자 하나당 효과를 추정하는 전 과정을 차례로 살펴본다. 각 방법이 데이터를 어떻게 해석하며, 행렬식을 이용해 $\beta$값을 계산하는지를 학습한다.
1. 일반 선형모형(GLM)의 행렬식
일반 선형모형(GLM)은 여러 요인(X)이 관찰값(Y)에 미치는 영향을 설명하기 위한 대표적인 통계모형이다. 행렬식을 표현하면 다음과 같다.
$$\mathbf{Y}=\mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon}$$
$\mathbf{Y}(n \times 1)$: 관찰된 표현형 벡터 (예: 개체별 체중)
$\mathbf{X}(n\times p)$: 디자인 행렬(어떤 개체가 어떤 요인에 속하는지 정리한 표)
$\boldsymbol{\beta}(p\times 1)$: 고정효과 벡터(요인별 효과 크기)
$\boldsymbol{\varepsilon}(n\times 1)$: 오차 벡터,
이 식을 풀어 $\beta$를 추정하려면, ‘추정방정식(estimation equations)’을 사용한다.
$$\mathbf{X}^T \mathbf{X} \hat{\boldsymbol{\beta}}=\mathbf{X}^T\mathbf{Y} \Rightarrow \hat{\boldsymbol{\beta}}=(\mathbf{X}^T \mathbf{X})^{-1}\mathbf{X}^T \mathbf{Y}$$
이제 실제 데이터에 이러한 방법들을 적용해 봅니다. 여기서는 한우의 근내지방함량(D_IMF)을 연구 대상 형질로, 약 100개체의 D_IMF 값과 후보 SNP 하나를 사용한다. 분석에 앞서 몇 가지 전 처리 및 확인 과정을 거친다.
• 데이터 로드 및 확인: .txt 형태의 데이터를 불러온 후, SNP 유전자형 데이터에 결측치(NA)가 없는지 확인한다. 실제 여기서 사용하는 데이터에는 결측치가 없다. 총 100마리의 한우(Hanwoo)데이터를 이용할 것이고, 표현형은 한우의 근내지방함량(D_IMF), 그리고 SNPx라는 유전자형을 이용하여 분석을 진행한다. 유전자형은 AA, AG, 그리고 GG형이며, AA형은 50마리, AG형은 30마리, 그리고 나머지 GG형은 20마리로 구성되어 있다.
• 유전자형 코딩 통일: 가끔 SNPx의 이형접합은 데이터에 AG와 GA 두 가지 형태로 기록되어 있는 경우가 있다. 이러한 경우 이는 동일한 유전자형을 의미하므로 둘 다 AG로 통일 처리해야 한다. R에서는 팩터로 읽을 때 자동으로 한 종류로 들어올 수도 있지만, 일반적으로 문자열 변환을 통해 통일시킨다.
• 집단별 근내지방함량 평균: 유전자형별로 근내지방함량의 평균과 표준편차를 구해본다. 예를 들어 SNPx의 경우 (결측 제외):
– AA형: 표본 수 50개, 평균 D_IMF ≈ 3.43%
– AG형: 표본 수 30개, 평균 D_IMF ≈ 4.01%
– GG형: 표본 수 20개, 평균 D_IMF ≈ 4.49%
이제 실제 데이터에 위 방법을 적용해 보겠습니다. 이번에는 한우의 근내지방도(D_IMF)를 표현형으로 사용하고, 특정 SNPx 마커 한 개를 선택하여 약 100개체 데이터를 분석합니다. Hanwoo_IMF.txt에는 다음 컬럼이 포함되어 있다.
• ID: 개체 식별자
• D_IMF: 근내지방도 (percent)
• SNPx: 선택한 SNP 유전자형 (AA, AG, GG)
• G: 대립유전자 a 개수 (0, 1, 2)
우리는 SNPx유전자형에 따른 한우의 근내지방함량의 차이를 분산분석으로 검정하고, 상가적 모형에 따른 효과도 추정해볼 것이다.
| 구분 | ANOVA | Regression |
|---|---|---|
| 가정 | 그룹 간 평균 차이 검정 | 대립유전자 개수에 비례 효과 |
| p-value | 7.8 × 10−27 | 5.7 × 10−28 |
| 효과크기 | AG–AA ≈ 0.58%p, GG–AA ≈ 1.06%p | allele 1개당 +0.54%p |
| 해석 | AG vs AA, GG vs AA 그룹차이 | 각 allele의 선형 기여 |
SNP 효과의 통계적 유의성
분산분석(ANOVA) 결과에서 $F(2,97)=119.0$, $p \approx 7.8 \times 10^{-27}$라는 값은, “세 유전자형(AA, AG, GG) 간 내장지방도의 평균 차이가 우연히 관측될 확률이 $10^{-26}$미만”임을 의미한다. 일반적으로 p<0.001 수준에서 ‘매우 유의하다’고 표현하므로, 이 SNP 하나만으로도 집단 간 내장지방도의 차이를 거의 완벽하게 설명할 수 있다는 강력한 증거이다. 또한, 분산분석의 결정계수 계산식
$$R^2 = \dfrac{SS_{SNP}}{SS_{Total}}=\dfrac{17.69}{17.69+7.21} \approx 0.710$$
를 통해, SNPx가 전체 분산 중 71%를 설명한다는 것을 확인할 수 있었다. 즉, 유전형이 한우 근내지방함량 변동의 절반 이상을 좌우한다는 것을 의미한다.
회귀분석 결과에서도 $t=15.43$, $p \approx 5.7 \times 10^{-28}$로, “G 대립유전자 한 개가 $D_{IMF}$를 변화시키는 선형 기울기(β=0.5378)가 우연히 0이 될 확률 역시 거의 0에 가깝다”는 것을 보여주었다. 여기서 절편(3.4424)은 기준 유전자형(AA)의 평균 근내지방함량, 기울기(0.5378)는 G 대립유전자가 하나 더해질 때마다 근내지방함량이 평균적으로 0.5378%p씩 증가한다는 실질적 효과 크기(effect size)를 제시한다. 또한 β의 95% 신뢰구간이 대략 [0.469, 0.607] 정도로 좁게 형성되어 있어(가령), “기울기가 0이 아니라는” 결론에 대한 불확실성이 매우 작음을 직관적으로 확인할 수 있다.
분산분석 후 post‐hoc contrast 혹은 직접 계산으로 “AA와 GG 평균 차이의 절반”을 상가적 효과로 산출하면
$$\dfrac{\mu_{GG}-\mu_{AA}}{2} \approx \dfrac{(3.432+2 \cdot 0.529)-3.432}{2}=\dfrac{1.058}{2}=0.529$$
가 되어야 하지만, 실제 데이터에서는 이 값이 0.529와 약간 다른 경우를 흔히 관찰합니다. 그 이유는 다음과 같다.
1. 표본 크기와 그룹 간 불균형 AG 집단이 30마리, GG 집단이 20마리로 표본 수가 다르므로, 분산분석의 “평균”과 회귀모형의 “OLS 평균 추정”이 동일하지 않을 수 있다. 회귀분석은 전체 100마리 데이터를 대상으로, 잔차 제곱합을 최소화하는 기울기를 찾기 때문에 각 그룹에 가중치(weight)가 다르게 반영된다.
2. 이형접합(AG) 평균이 완벽히 중간값이 아님 상가적 유전모형에서는 AG 평균이 AA와 GG 평균의 중간에 위치한다고 가정하지만, 실제 데이터에서는 AG 그룹의 표현형 분포가 순수 선형에 딱 떨어지지 않을 수 있다. 만약 AG의 실제 평균이 (AA+GG)/2보다 크거나 작으면, 분산분석에서 직접 계산한 μGG–μAA 절반과 회귀β는 미세하게 어긋난다. 유전학적으로 이는 우성의 효과가 0이 아님을 의미한다.
3. 회귀모형의 최소제곱 해(OLS)의 특성 OLS 회귀는 “잔차제곱합(residual sum of squares)을 최소화”하는 β를 산출하기에, 그룹별 평균 대신 전체 관측치의 공분산(cov(D_IMF, allele_G_count))과 분산(var(allele_G_count)) 비율
$$\hat{\beta} = \dfrac{\Sigma(x_i – \bar{x})(y_i – \bar{y})}{\Sigma (x_i – \bar{x})^2}$$
에 기반한다. 이때, 각 유전형에 속한 x값(0,1,2)의 분포와 y값의 편차가 균일하지 않으면, “단순히 평균 차이의 절반”과는 미묘하게 다른 값이 도출된다.
결국, 분산분석의 상가적 효과는 “두 동형접합 평균 차이의 이론적 절반”인 반면, 회귀분석의 β는 “실제 표본 데이터 전체에서 최소자승 조건을 만족하는 최적 선형 기울기”이기 때문에, 데이터 분포가 완벽히 이상적(이형접합이 정확히 중간, 표본 균형)이지 않으면 두 값은 미세하게 달라질 수밖에 없다.
* 통계적 유의성: ANOVA와 회귀 모두 $p<10^{-16}$ 수준으로 매우 유의하여, SNPx가 한우 근내지방함량에 미치는 영향이 통계적으로 확실하다.
* 효과 크기 비교: 분산분석 R²=0.71, 회귀분석 R²≈56으로, 두 모형이 설명하는 분산 비율에 차이가 있지만 모두 높은 설명력을 보였다.
* 모형 간 β 차이: 표본 크기 불균형, 이형접합 평균 위치 편차, 그리고 OLS 최적화 특성 때문에 “μGG–μAA 절반”과 “OLS β”는 완벽히 일치하지 않는다.
* 실무적 시사점: 이러한 미세한 차이까지 이해함으로써, SNP 효과 추정 시 “어느 방식이 더 적절할지”를 연구 맥락에 맞게 선택하고, 필요하다면 혼합효과모형이나 가중회귀 등을 이용해 보다 정교하게 모형을 개선할 수 있다.
분산분석과 회귀분석의 관계에 대해서는 통계 교과서와 GLM 이론에서 잘 설명되어 있으며, 상가적 유전모형과 마커 효과 추정에 대한 이론은 Falconer 등의 유전학 교과서에 상세하게 설명한다. 또한 본 예제에 사용된 한우 등심 근내지방함량과 유전자형 데이터는 공개 교육용 자료로서, 기존 한우데이터를 이용하여 본 실습에 맞게 100마리의 데이터셋을 시뮬레이션하여 제공하였다. 결과적으로 해당 SNPx는 한우 근내지방함량에 미치는 효과가 매우 크며, 이는 하나의 유전자가 표현형을 설명하는 설명력이 매우 크다고 할 수 있다. 그러나, 실제 육종현장에서 이러한 마커는 현실적으로 존재하기 어렵다.
본 교재의 내용 정리를 위해 연구논문 및 통계교재의 정보를 참고하였으며, 예제 데이터는 충남대학교 가축유전체육종학 연구실에서 제공하였다. 또한 회귀분석과 ANOVA의 관계에 대한 설명은 통계교과서를 참고하였다.
본인의 Google 계정으로 구글시트를 복사
본인의 Google 계정으로 구글코랩을 복사