CNVのDiploidコピー数情報から、コピー数アレル頻度をEMで推定する

ここ数日、Copy numer variationデータのケースコントロール解析処理についての覚書を続けている。昨日、CNVジェノタイプは、2本の染色体のコピー数の和として観測されることが多いこと、その情報から集団のコピー数頻度をEMアルゴリズムで推定すること、その頻度情報をもとにHWE検定をすることについてメモした。
今日は、コピー数の和としてのジェノタイプからアレル頻度をEM推定することについて、若干詳しい記載をする。

今、片方のアレルのコピー数がn1,もう片方のアレルのコピー数がn2とすると、この個人のコピー数の和はn1+n2である。
逆に、コピー数の和がNであるような情報が得られたとき、この個人が持つ、２アレルは、{0,N},{1,N-1},...,{N,0}のいずれかである。
コピー数がnkであるようなアレルの頻度を $p(nk)$ と書き表すと、コピー数の和がNと観測されたときに、個人の２染色体コピー数がn1,n2である尤度は、 $\frac{p(n1)p(n2)}{\sum p(ni)p(nj)}$ ただし、 $i,j \le N$ のように算出される。
今、この $p(nk)$ について、適当な初期値から開始して、仮に得られた $p(nk)$ を基にして、 $\frac{p(n1)p(n2)}{\sum p(ni)p(nj)}$ の割合で、コピー数Nの個体のアレル内訳を分配する処理を繰り返すと、いわゆるハプロタイプ頻度のEM推定と同様に収束する。ハプロタイプ推定のときにも、ローカルマキシマムへの収束が、まれにではあるが問題となるとおり、CNVコピー数アレル頻度のEM推定もそのリスクを負うが、実用上は(おそらく)無難な推定値をもたらすものと思われる。特に、HWEがおおまかになりたっているような集団に観測したデータの場合は大丈夫そうである。
推定がうまく行かないのは、コピー数アレルのうち、中ぐらいのコピー数のアレルの頻度が極端に小さく、２峰(多峰)性を示す場合などである。HWDにあるときには、２峰(多峰)性となりやすいことから、HWDのときの推定に誤差が出やすい。アレル数が２のときは、確定的に決まるので、推定結果が正確であることが、図から読み取れる。

図のファイル