EM 알고리즘 haplotype frequency 추정

EM 알고리즘을 통한 haplotype frequency 추정

우선 EM 알고리즘의 정의부터 알아보자.

EM 알고리즘은 관측되지 않는 잠재변수(unobserved latent variables)에 의존하는 확률 모델에서 최대우도(maximum likelihood)나 최대사후확률(Maximum A Posteriori)을 갖는 매개변수를 찾는 반복적인 알고리즘이다.

[출처] EM 알고리즘(1/2)|작성자 lhm0812

그렇다면 왜 haplotype frequency 추정할 때 EM 알고리즘 쓸 수 있는가? 관측되지 않은 변수가 있기 때문이다. 모수를 추정하고자 하는데 모수 추정에 필요한 값이 관측할 수 없는 잠재변수를 포함하고 있을 때 쓸 수 있다.

아래와 같은 세팅을 생각해보자.

Locus 1 has alleles A; a
Locus 2 has alleles B; b

Genotype counts are observed at the two loci for N individuals in the population.

homozygous for both loci: nAABB; naaBB ; nAAbb; naabb
homozygous for one loci: nAABb; naaBb; nAaBB; nAabb
heterozygous for both loci: nAaBb

근데 문제는 여기서 nAaBb는 2개의 가능한 하플로타입 nAB/ab or nAb/aB 이 있다. (즉, 두 loci 모두 genotype이 heterozygous인 경우 haplotype을 결정할 수 없다는 사실을 알 수 있다.)

우리가 알고 싶은건 pAB; pAb; paB; pab. 바로 이 haplotype frequency들이다.

여기서 pAB를 예로 들어 EM 알고리즘을 적용하는 과정을 살펴보자. 우선 pAB의 MLE는 아래와 같다.

p^AB = nAB/2N (N=관측된 사람 수. 하플로타입은 개인당 2개이므로 2N)

또한 nAB = 2nAB/AB + nAB/Ab + naB/AB + nAB/ab

이 때 문제는 nAB/ab를 모른다. 이를 잠재변수라고 하고, 이는 추정하려는 모수들 (pAB; pAb; paB; pab)과 관련되어 있다. 그래서 모수가 있다면 구할 수 있는 값이다. EM 알고리즘에서는 이를 모수를 적절한 값으로 초기화한 후, 반복적으로 nAB/ab를 구한다.

EM 알고리즘에서는 추정하려는 모수 pAB; pAb; paB; pab 를 적절한 값(ex. 0.25)으로 초기화한다. 그리고 이를 p0AB 처럼 표기한다.

관측한 값

Y = (nAABB; naaBB; nAAbb; naabb; nAABb; naaBb; nAaBB; nAabb; nAaBb)

와 초기화된 모수 p가 주어졌을 때, nAB의 기댓값은 아래와 같이 구할 수 있다.

여기서 E[nAB/ab]는 무엇인가?

왜 이렇게 되냐하면, nAB/ab는 nAaBb는 nAB/ab 또는 nAb/aB 이므로, 두 하플로타입(AB/ab, Ab/aB)이 등장하는 확률 값의 비율을 통해 nAB/ab의 기댓값을 구할 수 있다.

이를 통해 n0AB를 구했으면 pAB의 MLE를 통해 이를 업데이트한다.

p^AB = nAB/2N

이 과정을 반복하면 p^AB가 수렴해가고 이를 통해 추정값을 얻을 수 있다. 다른 모수들에 대해서도 이와 같이 한다.

출처 - http://courses.washington.edu/b516/lectures_2010/EM_Algorithm_Haplotype_Frequency_2010.pdf