Tuesday, March 3, 2015

Ancestry Composition in 23andMe

23andMe에서 Ancestry Composition 이라는 자신의 혈통적 뿌리가 어디서 유래 했는가를 알수 있는 서비스를 제공한다. 이는 아마도 23andMe 라는 회사가 기반을 둔 미국의 다인종 혼혈 상황에서 자신의 뿌리에 대한 사람들의 관심에 대한 결과가 아닐까 한다.

0. Contents

1. Basics

2. Wrinkle #1:  사람들은 일반적으로 여러 인종 혹은 혈통(popuy)을 조상으로 갖는다 (특히나 미국의 경우).

3. Wrinkle #2: 어떤 DNA 혹은 chromosome이 엄마, 아빠에서 온건지 알수 없다.

4. Ancestry Composition Overview

5. Prep 1: The DataSet

6. Prep 2: Population selection

7. Step 1: Phasing

8. Step 2: Window Classification

9. Step 3: Smoothing

10. Step 4: Calibration

11. Step 5: Aggregation & Reporting

12. Using Close Family Members

13. Testing & Validation

14. Ancestry Composition's Future




Basics

특정 DNA marker와 지리학적 위치의 연관성의 정보를 가지고 조상을 추측하는 것.
예를 들면 mtDNA haplogroup 의 한 종류인 haplogroup H 는 유럽에서 많이 발견되고 호주나 미국에서는 거의 발견되지 않는다. 만약 누군가 haplogroup H에 관련된 DNA marker를 가지고 있다면 그 사람은 유럽인일 확률이 크고 native american 확률은 작게 된다. 이렇듯 여러 DNA marker를 이용하여 혈통을 추측하는 것이다.

Wrinkle #1:  사람들은 일반적으로 여러 인종 혹은 혈통을 조상으로 갖는다

문제는 어떤이는 조상이 여러 위치에서 온 혹은 여러 인종이 뒤섞여 있을 수 있다는 것.
곧 admixture (the genetic mixing together of previously-separate populations), 원래는 유전적으로 분리가 되어 있던 집단이 서로 섞이는 것. 유럽인의 경우 조상이 유럽전역에 걸쳐 뒤섞인 경우.

Wrinkle #2: 어떤 DNA 혹은 chromosome이 엄마, 아빠에서 온건지 알수 없다.

chip 을 통한 genotyping은 unphased 상태, 곧 서로 다른 locus의 genotype이 같은 chromosome에서 왔지는 알 수 없는 상태이다. phased 상태의 SNP이 정보가 더 많다.

Ancestry Composition Overview

admixture and unknown phase 의 genotype 정보에서 혈통을 추측하기 위해

  1. 조상을 찾고자 하는  snp chip을 phasing을 통해 haplotype을 추정한다.
  2. 다음에 각 chromosome 을 특정 window 사이즈별로 겹치지 않게 구획을 나눠 reference genome과 비교하여 각 window 별 조상을 찾는다. 각 window는 여러 조상의 genotype이 섞인 것이 아니라 하나의 조상으로 부터 유래한다는 가정하이다. 
  3. "smoothing"을 통한 phasing correction
  4. 결과의 정확성을 위한 calibrate result 

Prep 1: The DataSet

10,000여 명의 데이터

  • 23andMe 자체 데이터 : 조부모, 부모의 출생지가 같은 회원의 정보
  • public resource : Human Genome Diversity Project, HapMap, 1000 Genomes project
reference 중 너무 genetic relationship이 가까운 개체는 제외. outlier로 판단되는 개체도 제외. 대략 10% 정보 제거됨


Prep 2: Population selection

위는 유럽인들의 유전정보만을 가지고 PCA를 그린 결과이다. 왼쪽 상단의 Finns 들의 경우 명확하게 isolation되는 것을 확인할 수 있다. 그러나 역사적으로 집단이 섞인 경우가 있어서 구분이 뚜렷하지 않은 경우도 있다. 이렇듯 구분이 되는 reference population 을 선택하였다.

Step1: Phasing

BEAGLE 프로그램을 customized 해서 Finch라는 프로그램을 개발, 이 프로그램을 통해서 haplotype phasing을 한다. BEAGLE은 phasing을 할 모든 데이터가 프로그램을 돌릴 때 갖춰져야 한다는 가정하에 있는데 추가적으로 고객이 발생하는 상황에서는 이 가정이 여의치 않기 때문에 이를 효과적으로 해결하기 위해 Finch를 개발.

Step2: Window Classification

phasing을 마친 후에 chromosome 을 100개의 marker 사이즈의 window로 쪼갠다. 보통 chromosome당 5,000 ~ 40,000 개의 marker가 존재. marker DNA 서열을 input으로 reference population에 의해 training 된 SVM 을 이용하여 각 window 의 ancestry를 결정한다.

Step3: Smoothing

SVM을 통해 ancestry 가 결정된 window의 정보는 first draft 이다. 여기에는 2가지의 에러가 포함 될 수 있다.

  1. window에 ancestry를 할당할 때 비슷한 정도의 ancestry candidate가 있을 경우 잘못된 ancestry를 할당하는 에러
  2. phasing step에서 haplotyping의 에러

위 두가지의 오류를 HMM (Hidden Markov Model)을 통해 수정하는 단계를 smoothing이라고 한다.

예를 들어 X, Y, Z 3개의 집단이 있고 SVM을 통해 아래와 같이 각 chromosome 의 window 별 ancestry 가 할당되었다고 할 때

parent1 의 뒤에서 두번째 window가  X라고 되어 있는데 이는 Z와 유사한 정도의 ancestry확률을 갖고 있어서 SVM이 잘못된 판단을 한 것일 수 있다. 또한 phasing step에서의 에러로 parent1 의 첫 4개 윈도우 X-X-Y-X는 parent 2의 Z-Z-Z-Z와 바뀌었을 확률이 높다.

이를 Hidden Markov Model 을 이용해서 correction을 하면 아래와 같게 된다.
위의 예는 smoothing을 설명하기 위해 단순화 한것이고 사실은 아래와 같이 각 window에 대해 모든 집단의 확률을 구한다.
위 그림은 African-American customer의 chromosome2에 대한 smoothing의 결과이다. 아래쪽 패널을 보게 되면 왼쪽부터 핑크-그린-핑크 (african-american-african)가 나오는 걸 확인할 수 있다. 

Step4: Calibration

위 단계에서 나온 그래프가 옳은지 우리는 어떻게 알 수 있을까? 시뮬레이션을 통해 systematic bias가 있는지 알아보고 이 systematic bias를 보정하는 recalibration단계를 추가한다.

Step5: Aggregation & Reporting

threshold를 정해서 그 threshold를 넘는 ancestry의 사이즈 에 비례해서 최종적으로 ancestry를 report하게 된다. 예를 들면 위 그림에서 오른쪽 부위의 green segment은 비록 이 region은 다른 집단에서 유래할 확률( orange segment :east asia)이 있더라도 70%의 threshold를 넘겼기 때문에 너비에 비례한 만큼 (전체 chromosome에서 0.26 % 차지) native americas 의 확률에 더하게 된다. 좌측 중간의 blue segment의 경우 threshold를 넘긴 특정 european population이 없기 때문에 report를 하지 않게 된다. 이런 경우에는 각 집단의 확률의 합이 threshold를 넘길 경우 이 집단들을 포괄하는 group 인 Broadly Northern European 으로 report를 하게 된다.

Using Close Family Members

부모나 자식의 데이터가 있다면 Finch를 통한 phasing 이 매우 정확해진다. 이는 higher-resolution result로 연결된다. 자식보다는 한명의 부모의 데이터가, 한명의 부모의 데이터보다는 두명의 부모 데이터가 있는 것이 조금 더 정확해진다.

13. Testing & Validation

14. Ancestry Composition's Future





* 위 내용은
https://www.23andme.com/ancestry_composition_guide/
http://blog.23andme.com/wp-content/uploads/2012/11/20121027_ancestry_painting_methods_poster.pdf
를 참조함

Wednesday, February 25, 2015

Predicting human height by Victorian and genomic methods

요즘 한창 관심 있는 분야가 breeding, 그 이유는 genomics를 이용한 prediction이기 때문이다.

인간이 본능적으로 관심이 가는 것은 아마도 생존과 번식일 것인데 NGS의 발전으로 한참 hot한 분야가 personalized genomics 혹은 cancer genomics가 아닐까 싶다. 이 분야는 생존과 관련된 분야로 개체의 건강과 수명에 관심에 의한 것이다. 나머지 하나가 번식. 좀 더 나은 자손 혹은 내 유전자의 생존의 확률을 높이고자 하는 본능. breeding이야 말로 genomics 가 번식에 공헌할 수 있는 분야라는 생각에 흥미를 가지고 있다.

사실 breeding은 가축 특히나 생장 기간이 긴 소의 경우 역사가 오래된 학문인데 그 연구의 근간은 통계혹은 수학적 모델링과 함께 한다. 특히나 EBV(estimated breeding value) 를 구하기 위한 BLUP (best linear unbiased model) 은 선형식에 대한 깊은 이해야 필요하다. 그렇기에 linear modeling 에 대한 공부를 진행하고 있는데 저 밑딴에 깊숙히 공부하다보며 어느덧 가던 길을 잃어버리고 내가 지금 무엇을 하고 있나하고 방황을 하게 된다. 그런의미에서 기술의 저 깊은 바닥에서 올라와 한참 적용이 되고 있는 상황이나 최신의 동향을 훑어보는 것이 좋은 방법이라 생각한다. 그런의미에서 오랜만에 논문을 리뷰해본다.

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를 바탕으로 함

Wednesday, January 21, 2015

Statistical models in human genetics (lesson 1)

Contents 
  • genetic mapping
  • data for genetic study
  • 3 common question in genetic mapping
  • common designs for genetic studies



Genetic Mapping
형질의 유전패턴과 chromosomal region의 유전 패턴을 비교하여 "positional cloning" 하는 것
- positional cloning : 어떤 유전자인지를 아는 것과 상관없이 어디에 위치하는지를 찾는 것.


Data for a Genetic Study
  • pedigree : 관계가 알려진 개체 모임 (예: 가족)
  • genetic markers : 손쉽게 측정 가능한 genetic variant
  • phgenotype : mendelian phenotype 혹은 complex phenotype

3 Common Question in Genetic mapping
  • 어떠한 형질에 영향을 주는 유전자가 있느냐? Is a trait genetic? => Epidemiological study
      • 집단과 친척관계에서 형질의 분포를 조사함.
      • Crohn's disease 의 경우 전체 인구에서의 발병률과 twins of affected individuals 차이가 남 
    • 그 유전자는 어디에 있느냐? Where are those genes? => Linkage analysis
      • Crohn's 이 있는 가계도와 genetic marker를 비교해 보니 "D16S3136" marker가 연관됨을 확인. "D16S3136" marker 의 위치는 알려져 있음.
    • 그 유전자는 무엇이냐? What are those genes? => Association analysis
      • Crohn's 환자 중 NOD2의 유전자의 mutation 을 갖고 있는 환자가 많음.

    Common Designs for Genetic Studies
    • parametric linkage analysis
      • 특정 모델와 위치를 평가하는 방식(disease loci의 allele frequency, genotype 별 발병률)
      • Identify variants with relatively large contributions to disease risk 
    • allele-sharing analysis
    • association analysis
      • Identify genetic variants with relatively small individual contributions to disease risk 
      • 가장 단순한 방법으로는 case와 control의 allele frequency 비교
      • genome-wide하게 수백만 개의 marker를 가지고 살수도 있고 몇개의 candidate에 초점을 맞춰서 할수도 있음




    * 위 내용은 http://www.sph.umich.edu/csg/abecasis/class/666.01.pdf 를 바탕으로 함

    Tuesday, January 20, 2015

    ipython notebook to pdf

    1. ipython nbconvert sample.ipynb --to latex --post PDF
    trouble: 에러 발생 (원인 불명)

    2. ipython nbconvert sample.ipynb --to latex 후 pdflatex sample.tex
    trouble : 한글의 경우 utf-8 인대도 에러 발생

    3. brew install Caskroom/cask/wkhtmltopdf  깔고
    ipython nbconvert sample.ipynb --to html 후 wkhtmltopdf sample.html sample.pdf
    trouble :  plot이 page 넘어가는 부분에 있을 경우 plot이 쪼개짐

    Friday, January 16, 2015

    Model of Evolution

    아래 내용은 http://www.evolution-textbook.org/content/free/contents/ch28.html 을 요약한 내용임

    Mathematical theory can be used in two distinct ways

    진화 생물학을 이해하는데 도움을 주는 수학적 기술은 두가지 방향으로 이용된다.



    1. 다양한 진화 과정의 결과를 산출 및 예측하고 진화 과정에서의 핵심 인자를 찾는데 이용된다 (위 그림 A, 점은 gene을, 색깔은 allele을 뜻한다.).
    •  여러 유전자의 영향을 받는 complex trait 가 어떻게 진화할것인가? 혹은 random drift, gene flow, selection이 어떻게 집단을 변화시킬것인가? 
    •  수학적 기술은 구체적이고 정확한 양적 예측을 위해 이용되기보다도 위 질문들에 대한 대답을 찾기 위한 이해를 돕는데 이용되어 진다. 
    • 수학적 기술을 이용한 "toy models"을 통해 최소한 원리적으로라도 언어로 주장되어 지는 것들의 효과를 명확하게 나타낼 수 있다.

    2. 거꾸로 시간을 거슬러 집단 내에서 샘플링된 개체 혹은 유전자의 기원을 찾는데 이용된다 (위 그림 B).
    • 범죄현장의 샘플에서 그 샘플의 인종을 추정한다던지,  유전병의 원인이 되는 변이의 나이를 추정한다던지.

    Deterministic processes

    A population is described by the proportions of the genotypes it contains
    집단이란 여러 타입의 entity의 집합이므로 각 타입의 비율을 나타내는 것이 중요하다. entity는 개체가 될수도 있고 유전자가 될수도 있다. 유전자가 entity일 경우 allele 이 유전자의 다른 타입이 되겠다. 그리고 집단에 대한 모델을 구축할때 variable과 parameter를 구분하는게 중요하다.

    • variable은 allele frequency와 같은 것
    • parameter는 selection coefficient, recombination rates, mutation rate와 같은 집단이 어떻게 진화할지슬 결정하는 요인.


    Reproducing populations tend to grow exponentially
    Gradual evolution can be described by differential equations
    populations tend to evolve toward stable equilibria
    stability is determined by the leading eigenvalue
    stability anlysis extends to several variables

    Random processes

    Our understanding of random events is relatively recent
    A random process is described by its probability distribution
    Probability can be thought of in several ways
    Probabilities can be assigned to unique events
    The Increase of Independently reproducing individuals is a branching process
     The sum of many independent event follows  a normal distribution
    Random walks can be approximated by diffusion

    Wednesday, January 14, 2015

    SNP chip analysis (from raw data to phylogenetic tree)


    아래 과정은 snp chip에서 case-control GWAS study까지의 대략적인 코드이다.


    Parse raw file

    describe parsing script.
    In [1]:
    #parse raw file
    
    import csv
    
    handle = open("snp_chips.csv")
    reader = csv.reader(handle)
    replicate = 4
    start_column = -16 * replicate + replicate - 1
    samples = reader.next()[start_column : : replicate]
    samples = ['snp_id'] + [sample.split('.')[0] for sample in samples]
    print 'samples      ::', samples
    print 'sample count ::', len(samples) - 1
    
    
    str_output_ori = '{}\n'.format(','.join(samples))
    str_output_val = str_output_ori
    str_output_val_filtered = str_output_val
    str_output_ori = 'chr,pos,' + str_output_ori
    pure_row_count = 0
    for sp in reader:
        if sp[0][0] == '#':
            continue
        snps = [sp[1]] + sp[start_column :: replicate]
        str_snps = ','.join(snps)
        str_output_ori += '{0},{1},{2}\n'.format(int(sp[53][-2:]),
                                                 sp[56], str_snps)
        str_snps = str_snps.replace("AAAA","0")
        str_snps = str_snps.replace("AAAC","1")
        str_snps = str_snps.replace("AACC","2")
        str_snps = str_snps.replace("ACCC","3")
        str_snps = str_snps.replace("CCCC","4")
        
        str_snps = str_snps.replace("AAAG","1")
        str_snps = str_snps.replace("AAGG","2")
        str_snps = str_snps.replace("AGGG","3")
        str_snps = str_snps.replace("GGGG","4")
        
        str_snps = str_snps.replace("TTTT","0")
        str_snps = str_snps.replace("TTTC","1")
        str_snps = str_snps.replace("TTCC","2")
        str_snps = str_snps.replace("TCCC","3")
        str_snps = str_snps.replace("CCCC","4")
        
        str_snps = str_snps.replace("TTTG","1")
        str_snps = str_snps.replace("TTGG","2")
        str_snps = str_snps.replace("TGGG","3")
        str_snps = str_snps.replace("GGGG","4")
        
        str_snps = str_snps.replace("X","nan")
        str_snps = str_snps.replace("NC","nan")
        str_output_val += '{}\n'.format(str_snps)
        if 'nan' not in str_snps:
            pure_row_count += 1
            str_output_val_filtered += '{}\n'.format(str_snps)
    
    #print 'snps which all samples have value :: ', pure_row_count
    
    handle.close()
    
    handle_w = open('snp_chips.parsed.csv','w')
    handle_w.write(str_output_ori)
    handle_w.close()
    handle_w = open('snp_chips.parsed.value.csv','w')
    handle_w.write(str_output_val)
    handle_w.close()
    handle_w = open('snp_chips.parsed.value.filterd.csv','w')
    handle_w.write(str_output_val_filtered)
    handle_w.close()
    
    samples      :: ['snp_id', 'sample1', 'sample2', 'sample9', 'sample3', 'sample5', 'sample14', 'sample13', 'sample11', 'sample16', 'sample12', 'sample15', 'sample6', 'sample8', 'sample7', 'sample4', 'sample10']
    sample count :: 16
    
    

    SNP simple stats

    In [2]:
    #SNP simple stats
    
    import pandas as pd
    import matplotlib.pyplot as plt
    
    data = pd.read_csv('snp_chips.parsed.value.csv', index_col=0)
    ## hist duplicated probes
    duplicated_snps = []
    indexes = list(data.index)
    for snp in indexes:
        duplicated_count = indexes.count(snp)
        if duplicated_count > 1:
            duplicated_snps.append(snp)
    print "duplicate probes :: ", len(set(duplicated_snps))
    
    data["index"] = data.index
    data.drop_duplicates(cols='index', take_last=True, inplace=True) # remove
                                                            # duplicate row
    del data["index"]
    
    ## plot non-nan count of samples
    snp_count_per_sample = data.apply(lambda x: sum(x.notnull()))
    fig, ax = plt.subplots(figsize=(15,4))
    ax.plot(snp_count_per_sample, marker='o', ls='-', color='blue')
    ax.set_xticks(range(0,16))
    ax.set_xticklabels(snp_count_per_sample.index, rotation=-45)
    ax.set_xlim([-1,16])
    ax.grid(True)
    ax.set_title('# of SNPs with value per sample')
    ax.set_xlabel('sample names')
    ax.set_ylabel('# of SNPs with value')
    for i in range(len(snp_count_per_sample.index)):
        ax.text(i, snp_count_per_sample[i], snp_count_per_sample[i], 
                ha='center', va='bottom')
    
    duplicate probes ::  455
    
    
    In [3]:
    data = data.loc[data.apply(lambda x: sum(x.isnull()),axis=1) == 0,:]
    print 'filtered data ::', data.shape
    
    filtered data :: (2494, 16)
    
    

    PCA dimension reduction

    In [4]:
    #PCA dimension reduction
    
    from sklearn.decomposition import PCA
    import math
    import numpy as np
    
    
    pca = PCA(n_components=3)
    pca.fit(data.T)
    tc = pca.transform(data.T)
    
    colors = cm.rainbow(np.linspace(0, 1, len(tc)))
    fig, ax = plt.subplots(figsize=(10,10))
    for i, sample in enumerate(data.columns):
        ax.plot(tc[i, 0], tc[i, 1], marker='o', ls='*', color=colors[i])
        ax.text(tc[i, 0], tc[i, 1], sample, ha='left', va='bottom', 
                color=colors[i])
    ax.set_title('PCA reduction')
    ax.set_xlabel('principle component 1')
    ax.set_ylabel('principle component 2')
    
    Out[4]:
    <matplotlib.text.Text at 0x107666510>
    

    Phylogenetic tree

    1. Calculate genetic distance

    • Calculate Genetic distance (http://www.plosgenetics.org/article/info%3Adoi%2F10.1371%2Fjournal.pgen.1002451)
      : Genetic distance (D) between pair-wise combinations of individuals was calculated as PLINK where

      \[ D = 1−{(IBS2+0.5 \times IBS1) \over N} \]
      IBS2 and IBS1 are the number of loci that share either 2 or 1 alleles identical by state (IBS), respectively, and N is the number of loci tested.
    • Make distance matrix file (phylip input file)
      : The input format for distance data is straightforward. The first line of the input file contains the number of species. There follows species data, starting, as with all other programs, with a species name. The species name is ten characters long, and must be padded out with blanks if shorter. And every cell should be space delimited.
    In [7]:
    ## calculate genetic distance
    import math
    
    distance_file = 'distance_matrix.txt'
    
    def calculate_distance(s1, s2, allele_num):
        return 1 - (allele_num - (s1 - s2).apply(math.fabs)).sum() / (allele_num * len(s1))
    distance_matrix = data.apply(lambda x1: data.apply(
                            lambda x2: calculate_distance(x1, x2, 4)))
    str_output = '{0}\n'.format(len(distance_matrix))
    for sample in distance_matrix.index:
        sample_name = sample[:10] if len(sample) > 10 else sample + ' '*(10-len(sample))
        str_output += '{0} {1}\n'.format(sample_name,
            ' '.join(['{:.4f}'.format(v) for v in distance_matrix.loc[sample,]]))
    
    handle = open(distance_file, 'w')
    handle.write(str_output)
    handle.close()
    

    1. Convert distance matrix to newick format

    Use Phylip (check "Running the programs in background or under control of a command file" in http://evolution.genetics.washington.edu/phylip/doc/main.html).
    In [8]:
    import os
    from IPython.display import Image
    
    if os.path.isfile('outfile'):
        os.remove('outfile')
    if os.path.isfile('outtree'):
        os.remove('outtree')
    str_output = '{0}\nY\n'.format(distance_file)
    handle = open('fitch_input','w')
    handle.write(str_output)
    handle.close()
    os.system('fitch < fitch_input > fitch_result')
    
    Out[8]:
    0
    

    1. Draw tree

    use phylip drawgram, drawtree
    In [9]:
    ## draw tree with python package "ete2"
    
    from ete2 import Tree, TreeStyle
    
    handle = open('outtree')
    lines = handle.readlines()
    handle.close()
    newick = ''.join([line.strip() for line in lines])
    t = Tree(newick)
    t.show()
    #t.render('myTree.png', w=150, h = 150, units='mm')
    #Image(filename='myTree.png') 
    
    No module named MySQLdb
    
    
     MySQLdb module could not be loaded
    
    
    In [10]:
    ## draw tree with phylip "drawtree"
    ## phylip "drawtree" make ps format file. Converting file format is needed to render image.
    
    loc_fontfile = 'font1' # check the location of fontfile
    handle = open('drawtree_input', 'w')
    handle.write('outtree\n{0}\nY\n'.format(loc_fontfile)) 
    handle.close()
    commands = []
    commands.append('drawtree < drawtree_input > drawtree_result') # use phylip "drawtree"
    commands.append('mv plotfile drawtree.ps')
    commands.append('convert drawtree.ps drawtree.png') # use "convert"
    for command in commands:
        os.system(command)
    Image(filename='drawtree.png')
    
    Out[10]:
    In [11]:
    # draw tree with phylip "fitch" result
    
    handle = open('outfile')
    lines = handle.readlines()
    handle.close()
    start_index_str = 'Negative branch lengths not allowed\n'
    end_index_str = 'remember: this is an unrooted tree!\n'
    print '<Unrooted Tree>'
    print ''.join(lines[lines.index(start_index_str)+1: 
                                        lines.index(end_index_str)])
    
    <Unrooted Tree>
    
    
            +--sample16  
          +-7 
          ! +--sample2   
        +-9 
        ! ! +---sample3   
        ! +-3 
      +-2   !  +--sample4   
      ! !   +-13  
      ! !      +---sample5   
      ! ! 
      ! +----sample6   
      ! 
      !    +-----sample7   
      !    !  
      ! +-12     +---sample8   
      ! !  !  +-11  
      ! !  !  !  +----sample9   
      ! !  +-10  
      ! !     !  +----sample10  
      ! !     !  !  
      ! !     +-14   +---sample11  
      4-1        ! +-6 
      ! !        ! ! ! +---sample12  
      ! !        +-5 +-8 
      ! !          !   +----sample13  
      ! !          ! 
      ! !          +----sample14  
      ! ! 
      ! +----sample15  
      ! 
      +----sample1   
    
    
    
    
    

    Significant diff SNPs

    use plink (consider all of heterozygosity as 1) up to manhattan plot PLINK case-control association test result("assoc") fields summary
    CHR     Chromosome
    SNP     SNP ID
    BP      Physical position (base-pair)
    A1      Minor allele name (based on whole sample)
    F_A     Frequency of this allele in cases
    F_U     Frequency of this allele in controls
    A2      Major allele name
    CHISQ   Basic allelic test chi-square (1df)
    P       Asymptotic p-value for this test
    OR      Estimated odds ratio (for A1, i.e. A2 is reference)
    
    In [12]:
    # make PLINK input files
    group1 = ['sample1', 'sample2', 'sample3', 'sample4', 'sample5', 'sample6']
    group2 = ['sample7', 'sample8', 'sample9', 'sample10', 'sample11', 
              'sample12', 'sample13', 'sample14', 'sample15','sample16']
    
    ## make MAP file
    data_plink = pd.read_csv('snp_chips.parsed.csv', index_col=2)
    data_plink = data_plink.loc[data.index,:]
    data_plink['id'] = data_plink.index
    data_plink.drop_duplicates(cols='id', take_last=True, inplace=True)
    data_map = data_plink[['chr', 'id', 'pos']]
    data_map.to_csv('data.map', sep=' ', index=False, header=False)
    
    ## make PED file
    data_plink = data_plink.drop(['chr', 'pos', 'id'], axis=1)
    data_ped = data_plink.apply(lambda x: ' '.join('{0} {1}'.format(
                                            snp[0], snp[-1]) for snp in x))
    str_output = ''
    for sample in data_ped.index:
        str_output += '{0} {1} {2}\n'.format(sample, 1 if sample in group1
                                             else 2, data_ped[sample])
    handle = open('data.ped', 'w')
    handle.write(str_output)
    handle.close()
    
    # execute plink
    plink_command = 'plink --file data --assoc --no-fid --no-parents '\
                                    '--allow-no-sex --no-sex --map3 --adjust'
    os.system(plink_command)
    
    Out[12]:
    0
    
    In [13]:
    # draw Manhattan plot
    import matplotlib.cm as cm
    
    plink_map = pd.read_csv('data.map', header=None, sep=' ', index_col=1)
    plink_map.columns = ['chr', 'pos']
    handle = open('plink.assoc')
    lines = handle.readlines()
    handle.close()
    plink_pvalue = pd.DataFrame([{'id': line.strip().split()[1], 
                'pvalue':line.strip().split()[8]} for line\
                in lines[1:] if line.strip().split()[8] != 'NA' ])
    plink_pvalue.index = plink_pvalue['id']
    plink_map = plink_map.loc[plink_pvalue.index, :]
    result = pd.concat([plink_map, plink_pvalue], join='outer', axis=1)
    chr_x_axis = result.pivot_table(cols='chr', aggfunc='max')
    chr_x_axis = chr_x_axis.loc['pos',].cumsum()
    x_axis = result.apply(lambda x: (x['pos'] + chr_x_axis[x['chr']-1])\
                          if x['chr'] != 0else x['pos'], axis=1)
    y_axis = pd.Series([-math.log(float(x), 10) for x in result['pvalue']])
    
    fig, ax = plt.subplots(figsize=(15,4))
    
    colors = cm.rainbow(np.linspace(0, 1, len(chr_x_axis)))
    for i, chr in enumerate(chr_x_axis.index):
        rt = result.loc[result['chr'] == chr, :]
        x_axis = rt.apply(lambda x: (x['pos'] + chr_x_axis[x['chr']-1])\
                          if x['chr'] != 0else x['pos'], axis=1)
        y_axis = pd.Series([-math.log(float(x), 10) for x in rt['pvalue']])
        
        ax.plot(x_axis, y_axis, ls='*', marker='o',color=colors[i])
    ax.set_xticks([0] + [i for i in chr_x_axis.values])
    ax.set_xticklabels(['chr{}'.format(i) for i in chr_x_axis.index],
                       rotation=-45)
    ax.grid(True)
    ax.set_title('Manhattan plot')
    ax.set_xlabel('chromosome (bp)')
    ax.set_ylabel('pvalue (-log10)')
    
    Out[13]:
    <matplotlib.text.Text at 0x10a23cdd0>
    

    SNPs diversity

    ROD

    \[ ROD = 1 - \pi_{cul} / \pi_{wild} \] \[ \pi = \Sigma_{ij}x_ix_j\pi_{ij} = 2* \Sigma_{i=1}^n\Sigma_{j=1}^{i-1}x_ix_j\pi_{ij} \]
    where \(x_i\) and \(x_j\) are the respective frequencies of the ith and jth sequences, \(\pi_{ij}\) is the number of nucleotide differences per nucleotide site between the ith and jth sequences, and n is the number of sequences in the sample.
    - calculation example : http://www.vcu.edu/csbc/lfsc520/week7/Week7LectureNotes.pdf

    selective sweep analysis\({}^1\)

    We used allele counts at SNP positions to identify signatures of selection in sliding 40-kb windows, for pools of sequence data. For each pool and SNP, we determined the numbers of reads corresponding to the most and least abundant allele (\(n_{MAJ}\) and \(n_{MIN}\)). For each window in each breed pool, we calculated a pooled heterozygosity score: \(H_p = 2\Sigma n_{MAJ}\Sigma n_{MIN}/ (\Sigma n_{MAJ} + \Sigma n_{MIN}){}^2\), where SnMAJ and SnMIN are the sums of \(n_{MAJ}\) and, respectively, \(n_{MIN}\) for all SNPs in the window. Individual \(H_p\) values were then Z-transformed as follows: \(ZH_p = (H_p - \mu H_p)/ \sigma H_p\). \({}_{1. http://www.nature.com/nature/journal/v464/n7288/abs/nature08832.html}\)
    In []: