Tuesday, February 17, 2015

Quantitative Genetics of Human Traits

0. Contents

1. Quantitative Genetics of Human Traits

2. The First Human Genetic Markers : ABO Blood Groups

3. Human Genomic Variation

4. Population Subdivision, Haplotype inference and linkage disequilibrium

5. Human Population Genetics

7. DNA Fingerprinting

8. Parentage Analysis

10. DNA Substitution Models

14. Linkage Analysis




1. Quantitative Genetics of Human Traits

1.1 Anthropometric variation

프랑스의 Alphonse Bertilon 등이 1800년대 후반에 anthropometrics 라는 것을 개발한다. 이는 요즘 말하는 continuous trait 인 physical variation 을 측정하기 위한 방법들이다. Francis Galton 은 부모에서 자식에게로 형질 유전의 일반적인 규칙을 찾기위해 이러한 continuous trait 의 variation을 분석하였다. 이후 이같은 연구는 R.A. Fisher에 의해 멘델의 법칙이 재조명 되면서 continuous trait의 유전에 대한 패턴을 설명하기까지 별 관심을 받지 못한다.  Fisher의 이론은 20세기 초반에 human disease genetics, animal breeding 에 핵심적인 역할을 하는데,  집단 유전학에서 형질에 영향을 주는 유전자에 대한 이해 없이 친척간의 trait variation의 패턴을 예측하는 것을 "quantitative genetics" 라 한다.

  •  continuous trait : discrete trait와 달리 수치화되는 형질
  • complex genetic trait : 여러 유전자의 영향을 받는 형질
  • 일반적으로 continuous trait 와 같이 수치화되는 형질은 그 원인이 여러 유전자의 영향에 의한것으로 complex genetic trait임 

1.2 Fisher's model

\(L\)개의 유전자들에 영향을 받는 특정 continuous trait 가 m 값으로 측정되었을 때, 각 유전자가 D 혹은 d의 allele 을 갖고 D를 갖는다면 \(+a/2\), d 를 갖는다면 \(-a/2\) 만큼 m 값에 영향을 미친다면
\[m=\Sigma_{i=1}^L x_i + e\]
\[ x_i = \left\{\begin{array}{r1} +a & if \; DD \\ 0 & if \; Dd \\ -a & if \; dd \end{array} \right. \]
\(e\)는 평균 0 이고 \(\sigma_{\epsilon}^2\) 분산을 갖는 환경에 의한 random effect
집단에서 이 trait의 평균과 genetic 분산을 구하면
\[ \bar{m} = \Sigma_{i=1}^L (aP_i - aR_i) = a\Sigma_{i=1}^L (P_i - R_i) \]
\[\sigma_G^2 = a^2\Sigma_{i=1}^L[P_i(1 - P_i) + R_i(1 + _2P_i - R_i)]\]
\(P_i, Q_i, R_i\) 는 i번째 gene 의 \(DD, Dd, dd\)의 allele frequency
* 위 분산에 대한 유도식 필요
total population variance  는
\[ \sigma_T^2 = \sigma_G^2 + \sigma_{\epsilon}^2 \]
이고 heritability of trait는 \[ H =  \frac{\sigma_G^2}{\sigma_G^2 + \sigma_{\epsilon}^2}\]
 trait 에 관련된 locus가 많을 수록 형질의 측정치 m의 분포는 정규분포에 가까워 진다.

1.3 Trait correlations between relatives

친척간의 trait 의 측정치를 이용해서 heritability 를 예측할 수 있다.
변수 쌍 (\(x_i, y_i\)) 의 correlation coefficient는
\[ \rho = \frac{[\frac{1}{n}\Sigma_{i=1}^nx_iy_i - \bar x \bar y](\frac{n}{n-1})}{\sqrt{\sigma_x^2\sigma_y^2}} \]
\[\bar x = \frac{1}{n}\Sigma_{i=1}^nx_i, \quad \bar y = \frac{1}{n}\Sigma_{i=1}^ny_i\]
\[\sigma_x^2= \frac{1}{n-1}\Sigma_{i=1}^n(x_i - \bar x), \quad \sigma_y^2= \frac{1}{n-1}\Sigma_{i=1}^n(y_i - \bar y) \]
Fisher 에 의해 부모 자식간의 correlation coefficient는 \[\rho_{FS} = \frac{1}{2}(\frac{\sigma_G^2}{\sigma_G^2+\sigma_{\epsilon}^2}) = \frac{1}{2}H\]
이므로 부모와 자식의 correlation  coefficient 를 알고 있다면 heritability \((H)\)를 예측할 수 있다.

1.4 Galton's mistake: regression to the mean

Galton이  자식들의 평균키와 부모의 평균키를 비교하는 내용을 1886년에 Nature지에 출판했다. Galton은 이 분석에서 부모의 평균키가 집단의 평균보다 높은 경우 자식의 키는 부모의 평균키보다는 작은 키가 되고 반대로 부모의 평균키가 집단의 평균키보다 작은 경우 자식의 키가 부모의 평균키보다는 크게 나온다는 사실을 발견했다. Galton 은 이를 부모 윗대의 조상이 보다 평균에 가까운 키를 갖기 때문에 어떠한 유전적 효과에 의해 이러한 현상이 나타난다고 생각하였다. 그러나 이는 단순히 통계적인 현상일 뿐이다.
키가 정규분포라고 이고 유전이 안된다고 가정하면 먼저 부모키를 랜덤하게 뽑고 자식 키를 랜덤하게 뽑았을 때 부모가 극단적으로 작을때 자식의 키도 극단적으로 작을 확률이 적어진다.

2. The First Human Genetic Markers : ABO Blood Groups

2.1 Introduction

1900년에 Karl Landsteiner 가 ABO가 혈액 타입을 발견한다(이로 인해 1930년에 노벨상을 받음). 혈액 검사는 A 혹은 B 의 antigen을 넣음으로서 확인한다(B 형 antigen을 넣었을 때 혈액의 응고가 일어나면 혈액에는 B antigen이 없었던 것이므로 B형이 아니게 된다.  AB형의 경우 어떤 antigen을 넣더라도 응고는 일어나지 않는다).
O 형의 경우 O 형을 나타내는 allele \((i)\) 가 다른 혈액형의 allele \((I^A, I^B)\)에 대해 recessive해서 \(ii\) allele 쌍만이 O형을 나타내기 때문에 \(i\) allele frequency는 O형 phenotype frequency \((f_O)\) 로 부터 쉽게 추정이 가능하다(\(p_i = \sqrt{f_i} = \sqrt{f_{ii}}\)). 그러나 phenotype frequency \(f_A, _B\) 로 부터 population allele frequency \(p_A, p_B\) 를 추정하는 건 쉽지 않다.
\[f_A = f_{AA}+f_{Ai}=p_A^2 + 2p_Ap_i =p_A^2 + 2p_A(1- p_A-p_B )\]
1955년 C.A.B Smith 에 의해 이 문제에 대해  expectation-maximization (EM) algorithm 이라고 불리는 방법을 이용한 해결방법이 고안되었다.



type O type A type B type AB
Observed 4,578 4,219 890 313

위와 같이 1000명의 사람들에서 혈액형을 조사했을 때 위 데이터로부터 EM algorithm 을 이요하여 allele frequency \(p_A, p_B, p_i \) 를 구해보도록 한다.

    1. 임의로 구하고자 하는 genotype의 개체수 (\(\hat n_{AA}, \hat n_{Ai},\hat n_{BB},\hat n_{Bi},\))를 정한다.
\[\hat n_{AA} = 100 \\  \hat n_{BB} = 100 \\ \hat n_{Ai} = n_{A} - \hat n_{AA} = 4219 - 100 = 4119 \\ \hat n_{Bi} = n_{B} - \hat n_{BB} = 890 - 100 = 790 \]
    2. "maximization" step : maximum likelihood 방식으로 현재 genotype count를 가지고 allele frequency를 계산한다.
\[ \hat p_A = \frac{2\hat n_{AA}+\hat n_{Ai}+ n_{AB}}{2n} = \frac{2 x 100 + 4119+313}{2 x 1000 } = 0.2316 \\ \hat p_B = \frac{2\hat n_{BB}+\hat n_{Bi}+ n_{AB}}{2n} = \frac{2 x 100 + 790+313}{2 x1000} = 0.06515 \\ \hat p_i = \frac{2n_{ii} + \hat n_{Ai} + \hat n_{Bi} }{2n}= \frac{2x4578+4119+790}{2 x 1000} = 0.70325 \]
    3. "expectation" step : Hardy-Weinberg Equilibrium 가정하에 2번 단계에서 추정한 allele frequency 를 이용하여 각 유전형의 수를 예측한다.
\[ \hat n_{AA}= n_A(\frac{\hat p_A^2}{\hat p_A^2+2\hat p_A\hat p_i}) = 596.4962  \\ \hat n_{Ai}= n_A(\frac{2\hat p_A \hat p_i}{\hat p_A^2+2\hat p_A\hat p_i}) = 3622.504 \\ \hat n_{BB}= n_B(\frac{\hat p_B^2}{\hat p_B^2+2\hat p_B\hat p_i}) = 39.40033 \\ \hat n_{Bi}= n_A(\frac{\hat p_B \hat p_i}{\hat p_B^2+2\hat p_B\hat p_i}) = 850.5997 \]
    4. 위의 "maximization" 과 "expectation" 단계를 반복한다.

2.2 Case-Control studies of disease-marker association

1925년 Bernstein 에 의해 tri-allelic locus for ABO blood types 가 발견된다. 1939년 Lionel Penrose 가 phenotypic trait간 linkage 를 검사하기 위해 association을 이용하는 일반적인 framework 를 만든다. 이후 1950년대에 초창기 case-control association study 중에 하나로 Aird 가 blood type과 5가지 암과의 association study를 진행한다 (아래 표, allele frequency 는 위 EM algorithm과 유사한 방법으로 구함).


2.2.1 Relative risk and odds ratio

Odds ratio (OR) 은 특정 인자가 기인하는 질병에 대한 위험도를 수치화하는 강력한 방법이다. 이는 Fisher (1935), Berkson (1953), Woolf (1955) 등 독립적으로 여러 연구에서 이용되었는데 Woolf 의 연구가 최초의 human genetic association analysis 의 적용이다.



case
(ex: cancer)
control
(ex: non-cancer)
exposed to the risk
(ex: smoker)
a b
unexposed to the risk
(ex: non-smoker)
c d

OR은 \[ OR =\frac{P_1}{1- P_1} / \frac{P_2}{1- P_2} = \frac{P_1(1-P_2)}{P_2(1-P_1)} \]
이고 종종 case와 control의 label 이 뒤바뀌더라도 값은 변하기 않고 기호만 바뀌기 때문에 Log Odds (LOD) 를 사용하기도 한다.
\[ LOD = log(\frac{P_1}{1-P_1}) - log(\frac{P_2}{1-P_2}) \]
\(P_1 = a/(a+c) =\) proportions of individuals exposed to the risk factor among cases
\(P_2 = b/(b+d) =\) proportions of individuals exposed to the risk factor among controls

*위 식에서 \(P_1, P_2\) 을 사용하면 보통 odds ratio 인 \(\frac{a}{b}/\frac{c}{d}\) 가 아니라 \(\frac{a}{c}/\frac{b}{d}\) 의 식이 된다. 하지만 결론적인 OR의 식은 동일하게 된다. 그리고 본문에 OR은 population frequency of the risk factor 에 독립적임을 증명하는 부분이 있는데 다시 한번 확인할 필요가 있다.

relative risk (RR)은 \[RR = \frac{Pr(case|exposed)}{Pr(case|unexposed)}=\frac{z_1}{z_2}\] 이고
\[OR=RR\times (\frac{1-z_2}{1-z_1})\] 이므로 \(z_1,z_2\) 가 매우 작으면(disease(=case)가 rare 한 경우) OR과 RR값은 거의 유사하게 된다.
보통 RR을 위의 식으로 바로 구할수 없는 경우는
\[RR = \frac{Pr(exposed|case)}{Pr(unexposed|case)} \times \frac{Pr(unexposed)}{Pr(exposed)} = (\frac{p_c}{1-p_c}) \times (\frac{1-p}{p})\]
\(p_c = a/(a+c)\) = frequency of the risk factor among cases
\(p = (a+b)/(a+b+c+d)\) = ovarall population frequency of the risk factor

2.2.2 Odds ratio estimators

population 전체에 대한 OR값을 추정하기 위해 sampling을 통한 OR 계산값을 사용한다. 이 OR값, 그러니까 관측된 frequency의 차이가 유의한가를 판단하기 위해 OR의 신뢰구간을 구해서 구해진 OR의 신뢰구간이 null hypothesis (OR = 1) 을 포함하고 있는지를 확인한다. 추정된 OR의 신뢰구간이 1을 포함하고 있으면 null hypothesis를 reject하지 못하고 1을 포함하고 있지 않으면 null hypothesis를 reject한다.
위의 표는 Aird (1954) 의 혈액형과 peptic ulcer 간의 linkage를 연구한 자료인다. 위 샘플에서 \(\hat{OR}\) 을 구하면
\[\hat{OR} = \frac{911\times 4219 }{579 \times 4578} = 1.45\]
이고 Woolf (1955) 가 고안한 OR의 standard deviation (\(\hat{\sigma}\)) 을 이용하면
\[\hat{\sigma} = \sqrt{1/a + 1/b + 1/c + 1/d} = \sqrt{1/911 + 1/579 + 1/4578 + 1/4219} = 0.057 \]
이고 95% 신뢰구간을 구하게 되면
\[\hat{OR} \in (1.45 - 1.96 \times 0.057 \times 1.45, \quad 1.45 + 1.906 \times 0.057 \times 1.45) = (1.28, 1.61) \]
로 이 95%신뢰구가는 1을 포함하지 않으므로 null hypothesis (혈액형과 peptic ulcer은 관련이 없다) 를 reject하게 된다.

3. Human Genomic Variation

3.1 Single Nucleotide Polymorphisms (SNPs)

mutation rate of nuclear DNA 는 굉장히 낮다(\(10^{-9}\)). 그래서 SNP이 있는 위치, 곧 locus에는 보통 2개의 allele(a set of variants)을 갖는다. 이 allele중 집단에서 비중이 적은 allele(=variant)를 minor allele 이라고 한다. 특정 locus (= location in DNA sequence) 에서 minor allele의 frequency가 predefined frequency(보통 0.05)보다 높으면 이 위치는 polymorphic 하다고 한다.

3.2 Genotypes, Haplotypes, Diplotypes

abc

3.3 Summarizing human genomic variations

abc

3.3.1 Allele frequency

abc

3.3.2 Genotype frequencies and Hardy-Weinberg proportions

abc

3.4 Genotype frequencies with inbreeding

abc

3.5 Population subdivision and the Wahlund effect

abc

4. Population Subdivision, Haplotype inference and linkage disequilibrium

4.1 Subdivision between Han chinese and Europeans
4.2 Marker heterozygosity
4.3 Inference of haplotype phase
4.3.1 Genotypes to diplotypes: A one to many mapping
4.4 Linkage disequilibrium

5. Human Population Genetics

5.1 Allele frequency change in populations
5.1.1 Continuous approximation of allele frequency
5.1.2 Migration between populations
5.1.3 Mutation
5.1.4 Selection
5.1.5 Genetic drift

7. DNA Fingerprinting

7.1 Introduction
7.2 DNA forensics
7.2.1 Calculating DNA fingerprint match probabilities

8. Parentage Analysis

8.1 Calculation of the likelihood ratio
8.2 Parentage probabilities using multiple loci

10. DNA Substitution Models

종간에 DNA substitution modeling에 대해 소개하고 species divergence에 대해 추정하고 sequence data로 부터 phylogenetic tree를 추측해본다.

10.1 Alignment of homologous sequences

phylogenetic analysis의 기본 전제는 sequence가 공통 조상으로 부터 유래하였다는 것이다. sequence alignment의 목표는 공통의 조상으로 부터 유래한 두개 혹은 그 이상은 sequence들의 유사부위를 정렬하는 것이다. 이 장에서는 alignment가 되어 있다는 가정하에 phylogeny를 추정하는 것에 초점을 맞춘다
* alignment는 Biological Sequence Analysis 라는 책을 참조한다


10.2 Pairwise percentage of substitutions

두 종의 molecular divergence 의 단순한 지표가 될수 있는 것은 두 sequence를 alingment한 뒤의
\[d = \frac{x}{n}\]
\(x\)는 두 sequence에서 차이가 나는 nucleotide의 갯수
\(n\)은 전체 비교 site

두 종간의 sequence 의 fixed difference를 DNA substitution 이라고 한다. 이 DNA substitution이 일어나는 이유는 한 종에서 mutation이 일어나고 이 mutation이 genetic drift 혹은 natural selection 에 의해 그 종에 fixed 되었기 때문이다.

mutation의 fixation에 영향을 미치는 것이 genetic drift 만 있을 경우(곧 neutral selection) Motoo Kimura 는 rate of substitution이 site-specific mutation rate와 동일함을 보였다. 이에 대한 증명은 아래와 같다.
diploid 인 chromosome을 갖는 개체의 특정 위치에서의 mutation rate 를 \(\mu\) 라 하면 size가 N (즉 개체수가 N) 인 집단에서의 예상되는 mutation의 수는 \(2N\mu\) 이다. 그리고 Sewall Wright에 의해 \(i\)개의 copy를 갖는 allele이 N개체 수의 집단에서 fixation 될 확률은 \(i/2N\) 임이 보여졌다.
rate of substitution, 곧 mutation 이 일어나고 fixation이 되는 rate는 (예상되는 mutation 갯수) x (fixation될 확률) 이라고 할수 있다 (mutation 이 fixation될 확률은 mutation은 하나의 새로운 allele이 나타나는 것이므로 \(1/2N\)이 된다).
\[v = 2N\mu \times \frac{1}{2N} = \mu\]
많은 genomic 영역에서의 mutation은 neutral하지 않은 경우가 많다. 이 경우 substution rate는 mutation이 neutral effect일 경우 보다 높거나 (positive selection) 낮게 된다(negative selection).
예를 들면 third codon position의 경우 이 위치에서의 nucleotide의 변화가 Amino Acid의 선택에 영향을 주지 않기때문에 이 위치에서의 substitution rate가 다른 위치의 substitution rate보다 높다. 이는 곧 Amino Acid의 변화가 negative selection임을 알 수 있다.

10.3 Modeling DNA substitutions

한 곳에서 substitution이 일어나고 같은 곳에서 또 substitution 이 일어날 경우 원래의 nucleotide로 돌아갈수 있기 때문에 단순히 substitution의 percentage 는 실제 일어난 substitution을 저평가 하게 된다. 이와 같은 문제를 처리 하기 위한 model 중 가장 초창기 모델이 Jukes and Cantor (1969), JC69 이다. JC69 model은 mutation 의 발생이 poisson분포임을 가정한다.
\[Pr(M) = \frac{e^{-vt}(vt)^{M}}{M!} \]
\(M\) = substitution 갯수
\(t\) = time
\(v\) = the rate of substitution per unit of time
* Poisson distribution :: \(f(k;\lambda)=Pr(X=k)= \frac{\lambda^ke^{(-\lambda)}}{k!} \)

위 식을 이용해서 substitution 이 한번도 안 일어날 확률은
\[Pr(M=0) = e^{-vt}\] 이고 한번 이상 substitution이 일어날 확률은
\[Pr(M \geq1)= 1-Pr(M=0)=1-e^{-vt} \]
또한 JC69 model은 4종류의 nucleotide로의 변화할 확률이 동일하다고 가정한다. 곧 특정 위치에서 substitution이 한번 이상일어나서 원래 T였던 allele이 A로 변화할 확률은
\[p_{TA}(t)=(1-e^{-vt})\frac{1}{4}\]
이다. 위 식의 앞부분은 한번 이상의 substitution이 일어날 확률이고 뒤의 \(\frac{1}{4}\)은 T allele이 A allele로 바뀔 확률이다.
반면바뀌거나 혹은 유지가 되서  T allele이 T allele로 될 확률은 아래와 같다.
\[ p_{TT}(t)=e^{-vt} + (1-e^{-vt})\frac{1}{4}=\frac{1}{4} + \frac{3}{4}e^{vt} \]

10.4 Substitution proportions under JC69 model

특정 위치에서 한번 이상의 substitution이 일어나서 최종적으로 원래의 allele과 다른 allele 될 확률은 \[ p_{i\not=j}(t) = (1-e^{-vt})\frac{3}{4}\]
이는 sequence의 모든 위치가 독립적이다라고 가정할 때 sequence에서 substitution이 일어나는 위치의 비율이라고 할 수 있다. 마찬가지로 sequence의 변화가 없는 비율은
\[p_{i=j}(t) = e^{-vt} + (1-e^{-vt})\frac{1}{4} = \frac{1}{4}+\frac{3}{4}e^{-vt} \] 가 된다.
위의 식들은 원래의 allele이 같은 allele로 변화 하는 것도 subtitution이라고 여겼다. 그러나 nucleotide가 변화가 없는 것은 substitution 이라 하지 않는다. 그렇기 때문에 원래의 substitution rate \(v\) 대신 그 값의 \(3/4\) 를 substitution rate라고 해야 JC69 model의 standard formulas가 된다.
\[\upsilon= \frac{3}{4}v\]이므로 \[v = \frac{4}{3}\upsilon\]를 위 식에 대체 해야 한다.

10.5 JC69 distance and divergence model

 nucleotide substitution의 비율 식을 변환하면 \(t\), 곧 시간에 관한 식이 되고 이는 두 sequence를 분리시킨 총 시간이 되고 두 sequence의 substitution rate 가 같다는 가정(molecular clock hypothesis)하에 이 total 시간은 divergence time의 2배가 된다는 사실로 divergence time을 추정할 수 있다.
위 nucleotide substitution proportion 식에서 \(v\) 대신 \(\frac{4}{3}\upsilon\)으로 대체한 뒤 \(t\)로 방정식을 풀게 되면 추정되는 total 시간은 아래와 같다.
\[\hat{t}=-\frac{1}{\upsilon}(\frac{3}{4})log(1-4/3\hat{p})\]
\(\hat{p}\) = the observed proportion of substitutions between sequences

substitution rate\(\upsilon\)을 안다면 \(t\)를 구할수 있고 divergence time 은 \(t/2\)가 된다.

예를 들어 사람과 침팬치의 non-coding nuclear DNA sequence의 substitution 비율이 대략 0.01~0.015가 되고 mammalian nuclear gene의 mutation rate \(\upsilon = 10^{-9}\) per year 정도로 알려져 있다. mutation이 neutral하다고 가정한다면 mutation rate가 곧 substitution rate가 되기 때문에 사람과 침팬치의 divergence time은 아래와 같이 계산 가능하다(\(\hat{p}=0.01 로 여김\)).
\[\frac{1}{2}\hat{t}=-\frac{1}{2} \times  \frac{1}{10^{-9}}(\frac{3}{4})log(1-4/3\times  0.01) = 5033633= 5 \mbox{MYA}\]

10.6 Kimura 2 parameter model

 JC69 model은 transition(A<->G, T<->C)과 transversion 이 일어날 확률을 동일하게 보기 때문에 현실적이지 않다. 여러 종의 sequence를 비교한 결과 transtion이 transversion 보다 훨씬 자주 일어난다. Kimura (1980)는 이에 대해 transition과 transversion의 rate를 다르게 하는 좀 더 현실적인 model을 제안했는데 이를 Kimura 2 parameter model (K80) 이라 한다.
첫번째 parameter가 relative rate of transitions versus transversions 인
\[\kappa = \frac{\alpha}{\beta}\]
\(\alpha = \) rate of transitions
\(\beta = \) rate of transversions

두번째 parameter인 ovarall substitution rate 는 아래와 같다.
\[\upsilon = \alpha + 2\beta\]
*transversion이 transition에 비해 2배의 가지수가 더 많기 때문에 x2를 한것으로 보임

위 두 parameter 하에 probability of transition
\[p_1(t) = \frac{1}{4} +\frac{1}{4}e^{-4\upsilon\kappa / (\kappa+2)} - \frac{1}{2}e^{-2\upsilon t(\kappa+1)(\kappa+2)}\]
이고 probability of transversion
\[p_2(t) = \frac{1}{4} - \frac{1}{4}e^{-4\upsilon\kappa / (\kappa+2)}\]
이다.
transition과 transversion의 expected proportion \(S,V\)가 \(S=p_1(t), V=2p_2(t)\) 이고
위 \(p_1(t), p_2(t)\)식을 \(S, V\) 치환한다음에 연립방정식을 풀게 되면
\[\hat{t}= -\frac{1}{\upsilon}(\frac{1}{2}log(1-2S-V) - \frac{1}{4}log(1-2V))\]
\[\hat{\kappa} = \frac{2llog(1-2S-V)}{log(1-2V)} -1\]
를 얻게 된다.

14. Linkage Analysis

linkage analysis의 기초 요소와 recombination rate, linkage mapping, disease gene mapping 을 추정하는 예를 살펴본다.

14.1 Probability model of recombination

 meiosis 단계에서 chromosome이 replicated 된 후에 chromosome segment를 교환해서 겹쳐진 상태인 chiasma를 형성하고 그 segment를 교환하는 recombination이 일어난다.
\(r_n\)을 n번 chiasma가 형성되었을 때 recombination이 일어날 확률이라고 하면 \(r_0 = 0\) 이고
\[r_n = \frac{1}{2}r_{n-1} + \frac{1}{2}(1-r_{n-1}) = \frac{1}{2}\]
중간식의 첫번째 부분 \(\frac{1}{2}r_{n-1}\) 은 n-1번째까지의 chiasma 형성 시 recombination 이 일어날 확률(\(r_{n-1}\))에 n번째 chiasma 형성에서는 recombination이 일어나지 않는(\(\frac{1}{2}\))  을 곱한 확률을 의미하고 두번 째 term인 \(\frac{1}{2}(1-r_{n-1})\) 은 n-1번째까지의 chiasma생성시 recombinationd 이 일어나지 않을 확률에 n번째 chiasma생성시에 recombination이 일어날 확률의 곱을 의미한다.

14.2 Haldane's map function


같은 chromosome  상에 있는 두 marker의 거리는 recombination 되는 생식세포의 비율(expected fraction of recombinant gametes, \(\theta\)) 로 측정될 수 있다.
\[\theta = \frac{1}{2}Pr(n>0)\]
\(Pr(n>0) =\) 1번 이상 chiasma를 생성할 확률 (위 식에서 \(\times 2\)를 하는 이유는 chiasma 형성된 뒤 교차가 일어날수도 있고 안일어날수도 있기 때문)로 Haldane 는 두 marker간의 거리는 멀고 chiasma의 생성이 일어날 경우가 매우 적기문에 이 확률은 Poisson 확률일 거라 가정
\[Pr(n) = \frac{e^{-cd}(cd)^n}{n!}\]
\(c = \) chiasma formation rate per Mb
\(d = \) marker A, B 의 Mb 거리
\[Pr(n > 0) = 1 - Pr(n=0) = 1 - e^{-cd} \]
이므로
\[\theta = \frac{1}{2}(1 - e^{-cd})\]
\(cd\) 가 작아질수록 \(\theta\)는 0 에 가까워 지고 이 상태의 두 marker를 complete linkage 라 한다. 반면 \( cd\)가 커질수록 \(\theta\)는 1/2에 가까워 지고 이 상태를 unlinked 라고 한다.

14.3 Inferring recombination rates

\(c\), chiasma formation rate 를 구하는 방법으로는 Taylor 급수를 이용한 linear appoximation이 있다.
Taylor 급수
\[e^{-x} = 1-x+\frac{x^2}{2} - \frac{x^3}{6} + \frac{x^4}{24} - ...\]
에서 \(x\)가 작으면 \(\frac{x^2}{2}\) 이 후의 식이 매우 작아지므로
\[e^{-x} \approx 1-x\]
가 된다. 고로
\[\theta =  \frac{1}{2}(1 - e^{-cd}) \approx \frac{1}{2}(1-[1-cd]) = \frac{1}{2}cd\]
위 식을 변형하여 chiasma formation rate 를 구하면
\[\hat{c}_1 = \frac{2\theta}{d}\]
가 된다. 이는 \(\theta \leq 0.01\) 일 경우에만 어느정도 추정값이 맞게 된다.
linear approximation말고 원래의 식으로 부터 바로 chiasma formation rate를 구하면
\[ \hat{c}_2 = \frac{1}{d}log(1-2\theta)\]
\(\theta\)의 단위를 cM (centi Morgan)을 많이 쓰는데 1 cM은 1% recombination per miosis를 의미한다. 그렇기 때문에 \(c\) 의 단위는 cM/Mb가 된다.

예를 들면,
500 그룹의 관계가 없는 부모와 자식한명 (unrelated family trio)의 genotyping을 하였을 때 두 marker의 총 500명의 자식의 genotyping이 되었기 때문에 가능한 haplotype은 1000개가 되고 이중에서 모계쪽에서 recombination 이 일어나서 자식에게서 발견된 recombinant haplotype (\(Y_m\)) 이 124,  부계쪽에서 recombination이 일어난 recombinant haplotype (\(Y_p\)) 이 86개 라고 하면
\[\theta = \frac{124 + 86}{1000} = 0.21 \]
두 마커의 거리가 1.6 Mb라고 할 때 linear approximation 에 의한 chiasma formation rate는
\[ \hat{c}_1 = \frac{2\times 0.21}{1.6} = 0.2625 = 26.25 \; \mbox{cM/Mb}\]
exact estimator는
\[  \hat{c}_2 = \frac{1}{1.6}log(1-2\times 0.21) = 0.340 = 34.0 \; \mbox{cM/Mb}\]
가 된다.

14.4 Linkage maps

abc

14.5 Linkage mapping of disease genes

linkage mapping  의 목적은 하나 이상의 genetic marker와 unobserved disease locus 간의 recombination proportion (\(\theta\)) 를 구하는 것이다.  Medelian genetic disorder는 하나 혹은 몇개의 genotype에 의해 일어나는 질병이기에 linkage mapping 이 효과적인다. disease locus 에 있는 genotype의 penetrance  는 \(Pr(disease|genotype)\) 으로 정의한다. \(D\)가 disease allele 이고 \(d\)가 normal alllele일때, 각 genotype의 penetrance는
\[f_1 = Pr(disease|DD) \\ f_2 = Pr(disease|Dd) \\ f_3 = Pr(disease|dd)\]
이고 simple recessive Mendelian disorder는 \(f_1 = 1, f_2 =0,f_3 =0\), 반면 simple dominant disorder \(f_1=f_2=1, f_3 =0\) 이다. \(f_3\)를 phenocopy rate라고 하는데 disease allele이 하나도 없을 경우의 발병률이다.

D allele 의 frequency 를 \(p_D\)라고 할때, HWE 가정하에
\[\begin{eqnarray}  f(DD) &=& p^2_D \\ f(Dd) &=& 2p_D(1-p_D) \\ f(Dd \times Dd) &=& f(Dd)^2 = (2p_D(1-p_D))^2 \\ f(Dd \times dd) &=& f(Dd) \times f(dd) = p_D(1-p_D) \times (1-p_D)^2 \end{eqnarray} \]
\(f(Dd \times Dd)\) 는 genotype \(Dd\) 를 갖는 두 사람이 mating을 할 확률이다. 매우 rare한 disease라면 \(p_D\) 값이 매우 작을 것이고 그러므로 \(f(Dd \times Dd)\), 곧 질병을 갖는 사람끼리 결혼할 확률과  \(f(DD)\), 곧 homozygous disease allele을 갖는 사람의 확률 또한 매우 작을 것이므로 고려 대상에서 제외해도 될 것이다. 곧 질병과 관련된 가계애서 한쪽 부모만이 heterozygous disease allele을 갖는 경우만 고려해본다.

unobserved disease locus (D:disease allele, d: normal allele) 와 genetic marker (A, a) 의 거리, 곧 recombination proportion 가 (\(\theta\)) 라고 하면 한쪽 부모의 genotype이 Aa/Dd, 다른 한쪽의 부모의 genotype이 aa/dd에서 생식세포의 종류와 확률은 아래의 표와 같다.

 그리고 질병을 가지고 있는 여러 가족의 자식의 genotype와 그 가계의 갯수가 아래 표와 같을 때
위와 같은 결과가 나올 확률은
\[\begin{eqnarray} Pr(\Upsilon_1, \Upsilon_2, \Upsilon_3, \Upsilon_4 \; | \; \theta) &=& \prod_{i=1}^4 Pr(\Upsilon_i \; | \; \theta) \\ &=& [\frac{1}{2}(1-\theta)]^{\Upsilon_1} \times [\frac{1}{2}\theta]^{\Upsilon_2} \times [\frac{1}{2}\theta]^{\Upsilon_3} \times [\frac{1}{2}(1-\theta)]^{\Upsilon_4} \\ &=& [\frac{1}{2}(1-\theta)]^{\Upsilon_1 + \Upsilon_4}\; [\frac{1}{2} \theta]^{\Upsilon_2 + \Upsilon_3} \end{eqnarray}\]
가 될 것이고 maximum likelihood 방식을 취하면 이 확률을 maximize하는 \(\theta\)가 disease locus 와 genetic marker의 거리가 될 것이다. \(\theta\) 값이 0이되면 두 loci는 complete linkage관계인 것이고 1\2가 되면 unlinked된 상태이다.
위 식의 log10을 위한 값을 lod score라고 한다.

A Primer on Linear Models (chapter 2)

Chapter 2 The Linear Least Squares Problem

2.0 contents

2.1 The Normal Equations
2.2 The Geometry of Least Squares
2.3 Reparameterization
2.4 Gram-Schmidt Orthonormalization
2.5 Summary of Important Results


2.1 The Normal Equations

One mathematical view of the linear model \(y = Xb + e\) is the best approximation \(Xb\) to the observed vector y.
Euclidean 방법으로 근사치를 구한다면 \[Q(b)= (y-Xb)^T(y-Xb) = \|y-Xb\|^2\] 를 최소화하는 \(b\) 를 \(Q(b)\) 의 least squares solution이라 한다.

2.2 The Geometry of Least Squares

2.3 Reparameterization
2.4 Gram-Schmidt Orthonormalization
2.5 Summary of Important Results

Monday, February 16, 2015

A Primer on Linear Models (chapter 1)

Chapter1 Examples of the General Linear Model

1.0 Contents

1.1 Introduction

1.2 One-Sample Problem

1.3 Simple Linear Regression

1.4 Multiple Regression

1.5 One-Way ANOVA

1.6 First Discussion

1.7 Two-Way Nested Model

1.8 Two-Way Crossed Model

1.9 Analysis of Covariance

1.10 Autoregression

1.11 Discussion


1.1 Introduction

\[y=Xb + e\]
\(y\) : N x 1 vector of observed reponses
\(X\): N x p matrix of fixed constants
\(b\): p x 1 vector of fixed, unknown parameters
\(e\): N x 1 vector of unobserved errors with zero mean

위 식이 general linear model 의 기본식. 1 장에서는 general linear model의 예를 알아본다.

1.2 One-Sample Problem

\[Xb = 1(\mu)\]
\(y_i\) 는 randomly sampled from a population with unknown mean \(\mu\) and variance \(\sigma^2\)

1.3 Simple Linear Regression

원래 식 \[y_i = \beta_0 +\beta_1x_i + e_i\]
general linear model 로 나타낸 식 \[y = \left[\begin{array}{c} y_1\\ y_2 \\ ... \\ y_N \end{array}\right], \quad Xb = \left[\begin{array}{cc} 1 & x_1 \\ 1 & x_2 \\ ... & ... \\ 1 & x_{N-1} \\ 1 & x_N \\ \end{array}     \right]\left[\begin{array}{c} \beta_0 \\ \beta_1 \end{array}\right] , \quad e = \left[\begin{array}{c} e_1\\ e_2 \\ ... \\ e_N \end{array}\right]  \]

1.4 Multiple Regression

원래 식 \[y_i = \beta_0 +\beta_1x_{i1} +\beta_2x_{i2} +\beta_kx_{ik} + e_i, \quad i = 1,...,N\]
general linear model 로 나타낸 식 \[y = \left[\begin{array}{c} y_1\\ y_2 \\ ... \\ y_N \end{array}\right], \quad Xb = \left[\begin{array}{cccc} 1 & x_{11} & x_{12} & ... & x_{1k} \\ 1 &  x_{21} & x_{22} & ... & x_{2k} \\ 1 & ... & ... & ... & ... \\  1 & x_{N1} & x_{N2} & ... & x_{Nk} \\ \end{array} \right] \left[\begin{array}{c} \beta_0 \\ \beta_1 \\ \beta_2 \\ ... \\ \beta_k \end{array}\right] , \quad e = \left[\begin{array}{c} e_1\\ e_2 \\ ... \\ e_N \end{array}\right]  \]

1.5 One-Way ANOVA

원래 식 \[y_{ij} = \mu + \alpha_i + e_{ij}\]
general linear model 로 나타낸 식 \[y = \left[\begin{array}{c} y_{11}\\ y_{12} \\ ... \\ y_{1n_1} \\ y_{21} \\ y_{22} \\ ... \\ y_{2n_2} \\ y_{a1} \\ y_{an_a}  \end{array}\right], \quad Xb = \left[\begin{array}{ccccc} 1_{n_1} & 1_{n_1} & 0 & ... & 0 \\ 1_{n_2} & 0 & 1_{n_2} & ... & 0 \\ ... & ... & ... & ... & ... \\  1_{n_a} & 0 & 0 & ... & 1_{n_a} \\ \end{array} \right] \left[\begin{array}{c} \mu \\ \alpha_1 \\ \alpha_2 \\ ... \\ \alpha_a \end{array}\right] , \quad e = \left[\begin{array}{c} e_{11}\\ e_{12} \\ ... \\ e_{1n_1} \\ e_{21} \\ e_{22} \\ ... \\ e_{2n_2} \\ ... \\ e_{a1} \\ ... \\ e_{an_a} \\ e_N \end{array}\right]  \]

\(i = 1,...,a\)
\(j = 1...,n_i\)

1.6 First Discussion

추후 정리 필요

1.7 Two-Way Nested Model

원래 식\[y_{ijk} = \mu + \alpha_i + \beta_{ij}+  e_{ijk}\]

\(i = 1,...,a\)
\(j = 1,...,b_i\)
\(k = 1,...,n_{ij}\)

\(\alpha \) 가 두 가지 F, P, 곧 \(\alpha_i = F, P\) 이고 각 \(\alpha\) 에 대해 \(\beta\) 가 \(D, M\) 과 \(A,B,E\),  곧 (\(\beta_F = {D, M}, \beta_P={A, B, E} \)) 일때


general linear model 로 나타낸 식  \[y = \left[\begin{array}{c} y_{F,D,1}\\ ... \\ y_{F,D,n_{11}} \\ y_{F,M,1} \\ ... \\ y_{F,M,n_{12}} \\ y_{P,A,1} /\\ ... \\ y_{P,A,n_{21}} \\ y_{P,B,1} \\ ... \\ y_{P,B,n_{22}} \\ y_{P,E,1} \\ ... \\ y_{P,E,n_{23}}  \end{array}\right], \quad Xb = \left[\begin{array}{cccccccc} 1_{n_{11}} & 1_{n_{11}} & 0 & 1_{n_{11}} & 0 & 0 & 0 & 0 \\ 1_{n_{12}} & 1_{n_{12}} & 0 & 0 & 1_{n_{12}} & 0 & 0 \\  1_{n_{21}} & 0 & 1_{n_{21}} & 0 & 0 &  1_{n_{21}} & 0 & 0 \\ 1_{n_{22}} & 0 & 1_{n_{22}} & 0 &  0 & 0 & 1_{n_{22}} & 0 \\ 1_{n_{23}} & 0 & 1_{n_{23}} & 0 & 0 & 0 & 0 & 1_{n_{23}}  \end{array} \right] \left[\begin{array}{c} \mu \\ \alpha_F \\ \alpha_P \\ \beta_{F, D} \\ \beta_{F, M} \\ \beta_{P, A} \\ \beta_{P, B} \\ \beta_{P, E}  \end{array}\right]  \]

1.8 Two-Way Crossed Model

원래 식 \[y_{ijk} = \mu + \alpha_i + \beta_j + \gamma_{ij} + e_{ijk}\]
\(i = 1,...,a\)
\(j = 1,...,b\)
\(k = 1,...,n_{ij}\)

위 식에서 \(\beta\) 가 없다면 이는 two-way nested model 이 되고 \(\alpha\)와 \(\beta\) 사이에 interaction이 없다면 \(\gamma\) 는 0이 된다.

\(\alpha\)와 \(\beta\) 사이에 interaction이 없고 \(n_{ij}\) = 1, (곧 no replication), \(i = 1,...,a\), \(j = 1,...,b\) 이며 balanced case 라고 한다면 맨 처음 식은
\[y_{ij} = \mu + \alpha_i + \beta_j +  e_{ij}\]

general linear model 로 나타낸 식은 \[y = \left[\begin{array}{c} y_{11}\\ y_{12} \\ ... \\ y_{1b} \\ y_{21} \\ y_{22} \\ ... \\ y_{2b} \\ ... \\ y_{a1} \\ y_{ab}  \end{array}\right], \quad Xb = \left[\begin{array}{cccccc} 1_b & 1_b & 0 & ... & 0 & I_b \\ 1_b & 0 & 1_b & ... & 0 & I_b \\ 1_b & 0 & 0 & ... & 0 & I_b \\ ... & ... & ... & ... & ... & ... \\  1_b & 0 & 0 & ... & 1_b & I_b \\  \end{array} \right] \left[\begin{array}{c} \mu \\ \alpha_1 \\ \alpha_2 \\ ... \\ \alpha_a \\ \beta_1 \\ ... \\ \beta_b \end{array}\right] \]

1.9 Analysis of Covariance

원래 식 \[y_{ij}\ = \mu + \alpha_i + \beta x_{ij} + e_{ij}\]
\(i=1,...,a\)
\(j=1,...,n_i\)

general linear model 로 나타낸 식은 \[y = \left[\begin{array}{c} y_{11}\\ y_{12} \\ ... \\ y_{1n_1} \\ y_{21} \\ y_{22} \\ ... \\ y_{2n_2} \\ ... \\ y_{a1} \\ ... \\ y_{an_a}  \end{array}\right], \quad Xb = \left[\begin{array}{cccccc} 1_{n_1} & 1_{n_1} & 0 & ... & 0 & x_1 \\ 1_{n_2} & 0 & 1_{n_2} & ... & 0 & x_2 \\ 1_{n_3} & 0 & 0 & ... & 0 & x_3 \\ ... & ... & ... & ... & ... & ... \\  1_{n_a} & 0 & 0 & ... & 1_{n_a} & x_a \\  \end{array} \right] \left[\begin{array}{c} \mu \\ \alpha_1 \\ \alpha_2 \\ ... \\ \alpha_a \\ \beta \end{array}\right] \]

만약 treatment \(\alpha\) 에 따라 covariate 의 slope 이 다르다면 식은 \[y_{ij}\ = \mu + \alpha_i + \beta_i x_{ij} + e_{ij}\] 가 된다.

* analysis of covariance는 식이 two-way nested model과 비슷하나 error가 \(e_{ijk}\) 가 아닌 \(e_{ij}\) 이다. 이는 오히려 하나의 factor 에 의해 샘플이 구분되는  one-way ANOVA 와 비슷한데 여기서 같은 그룹내에 각 샘플이 갖는 개별적인 속성(covariate) 의 변수를 넣은 것이라고 이해하면 될 것이다.

1.10 Autoregression

simple linear trend model for data collected over time \[y_i = \beta_0 + \beta_1t + e_t, \quad t = 1,2,...,N\]
\(i\)대신 \(t\)를 사용한 이유는 time ordering을 강조하기 위함.
추후 정리 필요


1.11 Discussion

위에서와 같이 다양한 통계 모델은 general linear model의 special case이다.  \(e_i\) 가 uncorrelated random variables with zero mean and constant variance \(\sigma^2\) 인 경우가 있고 다른 경우는 correlated 한 경우가 있다.
추후 정리 필요

Wednesday, February 11, 2015

PostgreSQL homebrew로 9.3 에서 9.4 로 upgrade 하기

Mac home-brew update & upgrade를 하면 PostgreSQL 이 버젼업 되는데 아래와 같은 에러가 난다.

LOG:  skipping missing configuration file "/usr/local/var/postgres/postgresql.auto.conf"
FATAL:  database files are incompatible with server
DETAIL:  The data directory was initialized by PostgreSQL version 9.3, which is not compatible with this version 9.4.0.
이때 아래 사이트를 따라하면 해결됨.

http://rcmachado.github.io/postgres/postgis/2014/12/20/upgrading-postgresql-using-homebrew.html

Wednesday, February 4, 2015

Linear Algebra

rank가 \(r\) 인 \( A_{m\mbox{x}n}\)  matrix 의 4가지 vector space 의 property 정리.

  • C(A) : column space of A. Linear combination of columns of A
  • R(A) : row space of A. Linear combination of rows of A
  • N(A) : null space of A. \(N(A) = \{y: Ay = 0 \} \)
  • N(A') : null space of A' (다르게 표현하면 \(\textrm{C}(A)\) 의 orthogonal complement). \(\textrm{C}(A) = S\) 일때 \(S^{\perp} = \textrm{N}(A')\) 


vector space C(A) R(A) N(A) N(A')
vector dimension \(R^m\) \(R^n\) \(R^n\) \(R^m\)
space dimension \(r\) \(r\) \(n-r\) \(m-r\)
projection matrix \(AA^-\) \(A'A'^-\) \(I - A^-A\) \(I - A'A'^-\)
- - - - -
- - - - -

* \(A'\) : transpose of A
* \(A^-\) : generalized inverse (GI) of a matrix A
* projection matrix : a square matrix \(P\) is a projection matrix that projects onto the vector space \(S\), \(S \subseteq R^n\)

  1. \(P\) is idempotene (=> \(PP = P\))
  2. \(Px \in S \; \forall x \in R^n \)
  3. \( Pz = z \; \forall z \in S \)



01PreliminaryLinearAlgebra

Fact V1
  • \(x_1,x_2,...x_n\) 이 LD (linearly dependent)  \( \Leftrightarrow \) \( x_j \) 는 \( x_1,x_2,..x_{j-1}\) 의 LC (linear combination)
Fact V2
  • vector space \(S\) 의 basis의 갯수가 n개이고 \(S\) 의 LI (linearly independent) 한  k개의 vector가 있다면 k \(\leq\) n
Fact V3
  • \({a_1,...,a_n}\) 과 \({b_1,...,b_k}\) 이 각각 vector space S의 basis라고 한다면 \(n=k\).
Fact V4
  • \(a_1,...,a_n\) 이 vector space S, S \(\in dim(\mbox{S})=n\) 의 LI vector라면 이는 곧 S의 basis 가 된다.
Fact V5
  • \(a_1,...,a_k\) 이 vector space S, S \(\in dim(\mbox{S})=n\) 의 LI vector라면 \(a_1,...,a_k\) 를 포함한 basis가 존재한다.
Fact V6
  • \(a_1,...,a_k\)이 \(R^n\)의 LI vectors 라면 \(a_1,...,a_n\) 이 orthornormal LI 가 되는  \(a_{k+1},...,a_n\) 이 존재한다.
Result A.1
  • \(rank(AB) \le min(rank(A),rank(B))\)
Result A.2
  • \(A=BC\) 이면 \(\textrm{C}(A) \subseteq \textrm{C}(A)\)
  • \(\textrm{C}(A)\subseteq \textrm{C}(B)  \) 이면 \(A=BC\) 를 만족하는 \(C\)가 존재한다.
Result A.3
  • \(rank(A_{mxn})=n \Longleftrightarrow \textrm{N}(A)={0}\)
TheoremA.1
  • \(rank(A_{mxn})=r\) 이면 \(dim(\textrm{N}(A))=n-r\) 이고 다르게 표현하면 \(dim(\textrm{N}(A_{mxn}))+dim(\textrm{C}(A_{mxn}))=n\)

02OrthogonalComplements

Result A.4
  • vector space \(S \subseteq R^n \) 이 있을 때 모든 \(x \in R^n\) 는 \(s\in S, t \in S^{\perp}\) 에 의해 \(x = s+ t\)  로 나타내어진다. 그리고 이 \(s,t\) 의 decomposition은  unique하다.
Result A.5
  • \(A_{mxn}\) 이 있을 때 \(\textrm{C}(A)^{\perp} = \textrm{N}(A')\)
Result A.6
  • \(S_1,S_2\) 는 \(R^n\) 안의 vector space 이고 \(S_1 \subset S_2\) 이면 \(S_2^{\perp} \subseteq S_1^{\perp}\) 이다.

04ProjectionMatrices

Result P.1
  • \(P\) 가 idempotent matrix (=> \(PP=P\)) 이고 vector space \(S\) 의 projection matrix 라면 \(S = \textrm{C}(P)\)   
Result A.14 
  • \(AA^-\) 는 \(\textrm{C}(A)\) 의 projection matrix
Result A.15
  • \(I - A^-A\) 는 \(\textrm{N}(A)\) 의 projection matrix 
Result A.16
  • \(P\) 가 symmetric 이고 idempotent라면 \(P\)는 \(\textrm{C}(P)\) 의 unique symmetric projection matrix 이다. 곧 \(P\) 가 자신의 column space 의 projection matrix가 된다.
Corollary A.4
  • \(P\) 가 symmetric idempotent matrix 일때 \(I - P\) 는 \(\textrm{C}(P)^{\perp} = \textrm{N}(P) = \textrm{N}(P') \) 의 sysmmetric projection matrix 이다.

05IntroGLM

asf

07PxNormalEquations

Theorem 2.1

  • \(P_x\) 는 \(\textrm{C}(X)\) 의 orthogonal projection 일 때, \(P_x = X(X'X)^-X'\) 이고 이는 아래와 같은 성질을 갖는다.
    1. \([X(X'X)^-X'][X(X'X)^-X'] = [X(X'X)^-X'] \)
    2. \( X(X'X)^-X'y \in \textrm{C}(X) \; \forall y \in R^n \)
    3. \( X(X'X)^-X'x = x \;  \forall x \in \textrm{C}(X)\) 
    4. \(X(X'X)^-X' = [X(X'X)^-X']' \)
    5. \(X(X'X)^-X'\) 는 \(X'X\)의 \( \forall GI\) 와 같다
Lemma 2.1
  • \(\textrm{N}(X'X) = \textrm{N}(X)\)
Result 2.1
  • \(\textrm{C}(X'X)=textrm{C}(X')\)
Result 2.3
  • \(Q(b)\)를 최소화 하는 \(\hat\beta\) 값이 \(\textrm{NE}(X'Xb=X'y)\) 의 해이다. 
Result 2.4
  • \(X'XA=X'XB \Longleftrightarrow XA=XB\)
Result 2.5
  • \((X'X)^-\) 가 \(X'X\)의 GI 이면 \((X'X)^-X'\) 는 \(X\)의 GI 이다. 곧 \(X(X'X)^-X'X = X\)
Corollary to result 2.5
  • \((X'X)^-\) 가 \(X'X\)의 GI 이면 \(X(X'X)^-\) 는 \(X'\)의 GI 이다. 곧 \(X'X(X'X)^-X' = X'\)
Result 2.6
  • \(I-P_x = I-X(X'X)^-X'\) 는 \(\textrm{C}m(X)^{\perp}=\textrm{N}(X')\) 
Theorem 2.2
  • \(\textrm{C}(W) \subseteq \textrm{C}(X)\) 일 때, \(P_x-P_w\)는 \(\textrm{C}((I-P_w)X)\) 의 orthorgonal projection 이다.
Normal Equations (NE)
  • \(X'X \hat\beta = X'y\) 
  • NE의 해는 Result 2.3 과 같이 \(Q(b)\)를 최소화 하는 \(\hat\beta\) 값이 되고 이는 \(X(X'X)^-X'y\) 이다.
  • solution of NE \( \Longleftrightarrow \hat\beta \Longleftrightarrow \mbox{b minimized of }Q(b) \Longleftrightarrow X(X'X)^-X'y\)
Corollary 2.1
  • NE 는 consistent 하다. 곧 항상 해가 있다.
Corollary 2.2
  • \(rank(X'X)=rank(X)\)
Corollary 2.3
  • NE의 그 어떤 해 \(\hat\beta\) 의 \(X\hat\beta\) 의 값은 동일하다.

08Reparameterization

\(y = W\alpha+\epsilon\) 과 \(y = X\beta+\epsilon\) 는 동등하다 혹은 reparameterization iff \(\textrm{C}(X)=textrm{C}(W)\)
Result 2.8
  • \(\textrm{C}(X)=\textrm{C}(W)\) 이면 \(P_x=P_x\).
Corollary 2.4
  • \(\textrm{C}(X)=\textrm{C}(W)\) 이면 \(\hat y = P_{x}y = P_wy \) 이고 \(\hat \epsilon = (I-P_x)y = (I-P_w)y = y- \hat y\)
Result 2.9
  • \(\textrm{C}(W)=\textrm{C}(X)\) 이고 \(\hat \alpha\) 가 NE \(W'W\alpha=W'y\) 의 해일 때, \(\hat \beta = T \hat \alpha\) 는 NE \(X'X\alpha=X'y\) 의 해이고  \(T \ni W= XT\)




Tuesday, January 27, 2015

Estimated Breeding Values (EBV)


EBV (Estimated Breeding Value) 란 상대적 유전값 (relative genetic value)으로 육종(breeding), 즉 인간이 원하는 방향으로 생물을 개량하는 데 객관적인 유전적 지표로 사용되어 진다.  아래 전체적인 EBV를 구하는 방법과 그 적용 예를 들어본다.

Contents
  • Introduction
  • Theoretical Background
  • Statistical Model
  • Computing and Using EBV's
  • How to use ebb's to make selection decisions
  • Multiple-trait model




Introduction

EBV 란 statistical prediction of the relative genetic merit, 곧 개체의 유전적인 유리함이 어떤지에 대한 통계적 예측값이다. 이 값은 육종에 사용된다. 원하는 형질의 EBV가 높은 개체만이 생식하게 나고 그렇지 않은 개체의 생식은 막음으로서 population 이 특정 형질에 특화되는 방향으로 흘러가게 하는 것이다.


Theoretical Background

표현형(phenotype) 은 고기의 지방함량, 사료대비 성장률등의 관측가능한 육체적 형질을 의미한다. 이 phenotype은 아래식과 같이 크게 유전적 요인과 비유전적요인으로 구성되어진다. 유전적 요인은 자손에게 유전이 되며 비유전적 요인은 환경적요인과 그밖에 그 개체에 존재하는 특성들을 의미한다. 비유전적 요인은 크게 fixed effect (sex, year of birth, management condition) 과 그밖에 측정하기 힘든 effect로 나뉜다.

\[ \mbox{Phenotypic observation = Environmental effects + Genetic effects + Residual effects} \]
이를 통계적 모델로 표현하면 아래와 같다.
\[ y_{ij} = \mu_i + g_i + e_{ij} \quad  (1)\]

\( y_{ij}\) = i번째 개체의 j 번째 관측치
\( \mu_i \) = i번째 개체의 fixed  environmental effect
\( g_i \) = i번째 개체의 유전형의 addictive(\(g_a\)), dominance(\(g_d\)), epistatic(\(g_e\)) genectic value의 총 합
\( e_{ij} \) = i 번째 개체의 j번째 관측에서의 random environmental effect의 총합

여기서 \(g_i\)를 (\(g_a\)), (\(g_d\)), (\(g_e\))로 대체하게 되면
\[\begin{eqnarray} y_{ij} & = & \mu_i + (g_a)_i + (g_d)_i + (g_e)_i + e_{ij}  \\  & = & \mu_i + (g_a)_i + e_{ij}^*    & (2) \end{eqnarray} \]
\(e_{ij}^* = (g_d)_i + (g_e)_i + e_{ij}\)

위 2번 식이 \(\mu_i\) 를 알고 있다는 가정하 EBV를 계산하는 기본식이다. 그러나  \(\mu_i\) 를 정확하게 알고 있는 것은 거의 불가능하기 때문에 같은 데이터를 가지고  \(\mu_i\) 와  \((g_a)_i\)를 동시에 구해야 한다. 이 두가지 값을 동시에 구하는 방법이 BLUP (Best Linear Unbiased Prediction) 이다.

BLUP은 동시에  fixed effect에 의해 발생하는 차이 조절하고 집단에서의 가계도에 있는 유전적 관계를 설명하는 과정에서 EBV를 구한다.
BLUP을 이용하기 위해서는 형질에 대한 heritability를 알고 있어야 한다.
BLUP을 표현하는 식으로 mixed model equation 이 있다.
\[ Y = X\beta + Zu + e \quad (3)\]
\(Y \) = vector of observed value
\(\beta \) = vector of fixed effect
\(u \) = vector of random genetic effect (2번식에서 \((q_a)i\) 와 같은 의미)
\(e \) = vector of residual error
\(X, Z \) = design matrix

3번식은 여러 단계를 거쳐 아래와 같은 4번식으로 유도가 된다.
\[ \left[\begin{array}{cc} X'X & X'Z \\ Z'X & Z'Z + A^{-1}\alpha \end{array}\right] \left[\begin{array}{cc} \beta \\ u \end{array}\right] = \left[\begin{array}{cc} X'Y \\ Z'Y  \end{array}\right] \quad (4) \]
\(A^{-1}\) = \(A\)의 inverse matrix,  A는 가계도를 바탕으로 개체간의 유전적 관계를 수학적으로 나타낸 matrix
\( \alpha \) = (1 - \(h^2\)) / \(h^2\), \(h\) 는 heritability
개체간의 유전적 관계에서 breeding value을 추측할 수 있는 BLUP의 핵심은 \(A^{-1}\) 에 있다.


Animal ID Sire ID Dam ID Sex Contemporary Group physical trait Test Score 
1 Unknown Unknown M
2 Unknown Unknown F
3 1 Unknown M
4 1 Unknown F
5 Unknown Unknown M
6 3 2 M 1 2
7 3 2 F 2 4
8 5 4 M 1 3
9 5 4 F 2 5

table 1

9마리의 개체가 있고 개체간의 가계도는 표1에서 와 같다고 하자. \(A\)는 개체간의 유전적 관계를 나타내는 matrix 로 9X9 diagonal matrix 이다.  \(A\)를 구하는 방법은 아래와 같다.

먼저 위의 표와 같이 가족 관계를 표를 만드는데 부모가 되는 개체가 자식이 되는 개체보다 앞이 되도록 표를 만든다.

for i to n,
    for j to i,

  1. 두 부모(s, d)를 모두 아는 경우 \[ \left\{\begin{array}{r1} if \: j < i,  a_{ij} = a_{ji} = 0.5(a_{js} + a_{jd}) \\ if  \: j = i, a_{ii} = 1 + 0.5a_{sd}  \end{array}\right. \]
  2. 한쪽 부모(s)만 아는 경우  \[ \left\{\begin{array}{r1}if \: j < i,  a_{ij} = a_{ji} = 0.5(a_{js} \\ if  \: j = i, a_{ii} = 1  \end{array}\right. \]
  3. 양쪽 부모를 다 모르는 경우 \[ \left\{\begin{array}{r1}if \: j < i,  a_{ij} = a_{ji} = 0 \\ if  \:   = i, a_{ii} = 1 \end{array}\right. \]


위 방법을 이용하게 되면 \(A\)는 아래와 같이 구해진다.
table2

위의 값이 클수록 두 개체간의 유전적 관계가 밀접한 것을 의미한다. 개체수가 많아질 경우 A row, column의 수는 기하급수적으로 많아지고 \(A^{-1} \) 의 계산은 더욱 힘들어 진다. 1970 중반에 \(A\) 를 구하지 않고 \(A^{-1} \) 를 바로 구할 수 있는 방법이 고안되면서 수십만 개체의 genetic relationship matrix 의 계산이 가능해졌다.


Statistical Model

BLUP이 fixed effect와 genetic effect를 동시에 고려하기 때문에 타당한 fixed effect 의 선택이 되어야만 개체간의 올바른 유전적 비교가 가능하다. fixed effect가 제대로 선택되지 않으면 fixed effect에 대한 효과가 residual error에 포함되어지고 그만큼 모델의 설명력은 떨어지게 된다.  fixed effect에는 sex, 그리고 physical interpretation을 수행한 평가자, 나이, 평가 시기 등이 있다.


Computing and Using EBV's

heritability(h) 를 20%, fixed effect 를 contemporary group이라는 가정하에 예제 표1의 EBV를 계산하면 아래와 같다.
\[ Y = X\beta + Zu + e \]

\[ \left[\begin{array}{c} y_6 \\ y_7 \\ y_8 \\ y_9  \end{array}\right] = \left[\begin{array}{c} 2 \\ 4 \\ 3 \\ 5  \end{array}\right] = \left[\begin{array}{cc} 1 & 0 \\ 0 & 1 \\ 1 & 0 \\  0 & 1 \end{array}\right]  \left[\begin{array}{c} b_1 \\ b_2   \end{array}\right] +  \left[\begin{array}{r1} 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0  \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0  \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0  \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1   \end{array}\right] \left[\begin{array}{c} a_1 \\ a_2 \\ a_3 \\ a_4 \\ a_5 \\ a_6 \\ a_7 \\ a_8 \\ a_9  \end{array}\right]  + \left[\begin{array}{c} e_6 \\ e_7 \\ e_8 \\ e_9  \end{array}\right]\]
가 된다.
\[ \left[\begin{array}{cc} X'X & X'Z \\ Z'X & Z'Z + A^{-1}\alpha \end{array}\right] \left[\begin{array}{cc} \beta \\ u \end{array}\right] = \left[\begin{array}{cc} X'Y \\ Z'Y  \end{array}\right] \quad (4) \]
을 이용하기 위해서 \(A^{-1}\)  와 \( \alpha \)를 구하게 되면

 \(A^{-1}\) 은 아래와 같고
table 3

heritability = 0.2 로 가정하였기 때문에 \( \alpha \) = (1 - 0.2) / 0.2 = 4 가 된다.

4번식에 \(X, Z, A^{-1}, \alpha \) 값을 넣고 계산하게 되면 아래와 같이 fixed effect의 값과 개체의 EBV를 구하게 된다.              
\[ \left[\begin{array}{c} \beta \\ u \end{array}\right] = \left[\begin{array}{c} \mbox{Contemporary group 1 average} \\ \mbox{Contemporary group 2 average} \\ EBV 1 \\ EBV 2 \\ EBV 3 \\ EBV 4 \\ EBV 5 \\ EBV 6 \\ EBV 7 \\ EBV 8 \\ EBV 9 \end{array}\right] = \left[\begin{array}{c} 2.5 \\ 4.5 \\ 0 \\ -0.09 \\ -0.07 \\ 0.07 \\ 0.09 \\ -0.13 \\ -0.13 \\ 0.13 \\ 0.13 \end{array}\right]  \]

How to use EBV's to Make Selection Decisions

EBV 는 그저 germplasm(= collection of genetic resources) 일 뿐 EBV가 높다고 항상 좋은 것은 아니다. phenotype을 disease susceptibility 로 두었다면 EBV 이 최소값을 같는 개체를 선발하는 것이 옳을 것이다.
위의 예의 physical trait를 사료 대비 성장률이라고 성장률을 높이길 원한다는 가정하면 먼저 성별로 개체를 나누고 EBV를 내림차순으로 정렬한다 (표 4).

Male ID Male EBV Female ID Female EBV
8 0.13 9 0.13
5 0.09 4 0.07
1 0.00 2 -0.09
3 -0.07 7
6 -0.13
table 4

trait 을 더욱 강화하고자 한다는 가정하에서는 male은 8번 개체를 female 개체는 9번 개체를 선택한다.



* 위 내응은 https://www.seeingeye.org/images/HowToUseEBV.pdf 를 바탕으로 함.

GWAS 2

일전에 GWAS (Genome-Wide Association Study)  에 대해 간략하게 포스팅을 한적이 있었는데 이번에 조금 더 자세하게 정리해 보려 한다.


Content
  • Important Questions in Human Genetics
  • Concepts Underlying the Study Design
  • Capturing Common Variation
  • Genotyping Technologies
  • Study Design
  • Association Test
  • Replication and Meta Analysis


Important Questions in Human Genetics

  • human genetics 의 핵심적인 목표는 common complex disease의 유전적 위험 요소를 찾는것. 여러 기술과 study design이 있는데 그중 하나가 GWAS. 
    • 유전적 위험 요소를 찾는 GWAS의 적용사례 중 하나가 age-related macular degeneration (AMD) 의 원인이 되는 Complement Factor H gene을 찾은 것
  • GWAS를 통한 personalized medicine 분야의 적용도 또다른 human genetics의 목표라 할 수 있음
    • GWAS를 통한 항응고제인 wafarin 복용량 결정이 한 예


Concetps Underlying the Study Design


  • 가장 일반적인 SNP을 이용함. SNP 은 일반적으로 MAF (minor allele frequency) 가 높은 single base substitution 을 뜻함
  • rare disorder 의 경우 특정 소수 혹은 단일의 유전자에 의한 것이 많기 때문에 몇몇의 가계도를 분석하는 linkage analysis를 통해 genetic risk factor 가 잘찾아짐. 하지만 common disease의 경우 linkage analysis 가 효과가 없었음
  • common disease는 다르게 접근해야 한다는 생각과 common disease의 원인 유전자가 몇몇 common variant (MAF가 높은 SNP) 이라는 발견과 맞물려 common disease의 원인은 common variant 라는  CD/CV (common disease / common variant) 가정이 성립됨. 이를 바탕으로 하는 것이 GWAS
    • CD/ CV 에 서 분파된 하나의 생각은 common variant가 disease에 영향을 주기 때문에 그 영향의 크기는 상대적으로(rare variant to rare disease에 비해) 적을 것이다(=low penetrance)
    • common variant가 low penentrance이지만 질병이 유전된기때문에 multiple SNP이 원인이 되는게 아닐까 생각하게 됨.
    • 위 두가지 생각으로 family-based에서 population-based study로 초점이 옮겨가게 됨
    • SNP의 allele frequency와 그 SNP에 의한 질병의 위험도(effect size)가 중요 요소가 genetic study design의 핵심. 위 그림에서 우측 하단이 GWAS 연구하기 좋은 케이스 
    • The National Human Genome Institute GWAS catalog 에 GWAS 를 통한 disease causal SNP 이 정리 되어 있음


Capturing Common Variation

CD / CV 가설을 제대로 테스트 하기 위해서는 몇가지 체계적 접근이 필요함

  • Common variant (SNP) 의 파악 필요. 
    • 연구해야 할 SNP들이 어떤것이 있는지 또한  genome 상의 위치 의 정보 필요 
  • population 간 SNP의 차이에 대한 이해 필요
    • population 간에 차이를 고려 하지 않으면 GWAS연구에서 target trait 관련 SNP이 아니라 population 간의 차이가 나는 SNP이 뽑힐 수 있음
  • SNP간의 correlation 정보 필요
    • 중복된 정보를 제거하기 위해 혹은 효율적인 연구를 위해
  • 위 3가지 목적에서 생긴 것이 The International HapMap Project
    • 1.6 M SNPs
    • 11 human population
    • LD (linkage disequilbrium) 측정
      • SNP 이 chromosome의 동일선상에 있을 경우 연관되어 있는 정도를 나타내는 것
      • african이 generation 이 다른 population에 비해 높기 때문에 LD 영역이 상대적으로 적음
      • recombination 이 일어나면 두 SNP은 LD가 깨짐, 곧 linkage equilibrium
      • LD의 지표로 \(D' \) 와 \( r^2\) 를 많이 사용
        • \(D'\)는 0 ~ 1 사이값을 갖음. 1 이 complete LD를 의미, 0 은 HWE 상태 하에 두 SNP이 독립일 경우.
        • \[ D' = \left\{\begin{array}{r1} \frac{\pi_{AB}\pi_{ab} - \pi_{Ab}\pi_{aB}}{\mbox{min}(\pi_A\pi_b,\pi_a\pi_B)}   & if \pi_{AB}\pi_{ab} - \pi_{Ab}\pi_{aB} > 0 \\ \frac{\pi_{AB}\pi_{ab} - \pi_{Ab}\pi_{aB}}{\mbox{min}(\pi_A\pi_B,\pi_a\pi_b)}   & if \pi_{AB}\pi_{ab} - \pi_{Ab}\pi_{aB} < 0   \end{array} \right\} \]
        • \( r^2\) 은 correlation 식
        • \[r^2 = \frac{(\pi_{AB}\pi_{ab} - \pi_{Ab}\pi_{aB})^2}{\pi_A\pi_b\pi_a\pi_B} \]
        • LD 분석을 통해 tagSNP 이 선택되어짐(LD block내 대표 SNP). population마다 LD가 다르기 때문에 tagSNP도 달라짐
        • 이 LD로 인해 GWAS로 뽑힌 결과가 꼭 질병과 direct하게 연관되어 있다는 보장은 없음. 뽑힌 SNP의 LD 관계에 있는 SNP이 질병의 직접적인 원인 SNP 일 수 있음.



Genotyping Technologies

Illumina 와 Affymerix 사의 SNP chip이 많이 사용됨. 요즘은 sequencing으로 많이 함.


Study Design

case control design 과 quantitative design

  • quantitative design
    • target trait (혹은 disease)의 지표가 quantitative 함. 곧 양적으로 측정됨
    • 통계적으로 유의한 genetic effect를 파악하기 쉬운 경향이 있음
    • heart disease의 지표가 될수 있는 HDL, LDL 등의 lipid level이 quantitative trait의 예
  • case control design
    • trait가 affected VS unaffected 식으로 이원화 된 경우의 study



Association Test 

사실상 GWAS분석이라는 것은 SNP 하나하나를 독립적으로 phenotype과 연관성이 있나 검사하는 것임

  • quantitative design
    • GLM(depedent variable를 함수로 변형시켜 linear 하게 만든 것) 또는 ANOVA로 분석
  • case control design
    • contingency table method(chi-square test, Fisher's exact test), logistic regression을 사용
    • logistic regression은 covariates 의 보정을 할수 있어 contingency table 보다  더 많이 쓰임

genotype data를 숫자로 변경하는 "data encoding"  하는 방식 선택도 테스트의 통계적 검정력에 영향이 있음(genotype-based group의 숫자에 따라 자유도 변화등).

  • allelic VS genotypic
    • allelic : SNP 의 하나의 allele과 phenotype의 연관성을 봄
    • genotypic : 한쌍의 allele 을 하나의 값으로 치환함. 
      • dominant, recessive, multiplicative, addictive model이 있음.
sex, age, study site, population substructure등의 covariate 보정필요

  • population substructure를 고려해야 하는 이유 : 예를 들어 폐암 환자그룹에 특정 population이 많이 있다면 폐암 관련 유전자가 아니라 population specific한 SNP이 뽑힐것임.
  • STRUCTURE나 ENGENSTART 방법으로 보정 필요. HapMap의 인종별 allele freuncy를 참고하거나 아니면 PCA를 이용.

SNP하나하나에 대한 통계적 분석이기 떄문에 multiple test correction 필요

  • Bonferroni correction : 모든  SNP이 independent 하다는 가정
  • FDR (false discover) 
  • permutation testing 
  • genome-wide significance : 특정 population 의 LD의 분포를 바탕으로 독립적인 "effective" genetic region 이 있음을 바탕으로  correction 을 해야하는 test 갯수를 기준으로 p-value threshold를 정함, 유럽인의 경우 7.2e-8


Replication and Meta Analysis


  • replication : GWAS 결과에서 얻어진 SNP을 검증하는 실험
  • meta analysis : 여러 study를 합쳐서 분석하는 방법



* 위 내용은 http://www.ncbi.nlm.nih.gov/pubmed/23300413 과 http://bioinformatics.org.au/ws09/presentations/ 의 day3 jstankovich의 pdf를 바탕으로 함