多人数一括ＤＮＡプロファイリングのための 確率計算法に関する考察

Presentation on theme: "多人数一括ＤＮＡプロファイリングのための 確率計算法に関する考察"— Presentation transcript:

1 多人数一括ＤＮＡプロファイリングのための 確率計算法に関する考察
法数学勉強会 2011/09/10 京都大学大学院医学研究科 統計遺伝学分野 奈良原舞子　山田　亮

2 状況 大災害が発生 多数の行方不明者 多数の身元不明遺体 外見や所持品などから身元が特定された遺体はすでに遺族に返還されている。
残っている遺体の手がかりは、主に遺伝情報

3 使えるデータ 遺体 ジェノタイプ 発見時の状況など 行方不明者 家系情報 家族のジェノタイプ

4 個人の鑑定 簡単に描くと 行方不明者　Missing 身元不明者　found Body m1 b1 集団の誰か

5 集団の鑑定 .. .. .. .. .. .. 行方不明者 Missing 身元不明者 found Body m1 b1 m2 b2 m3
mN bN

6 集団の鑑定 N人をN体に割り付ける N! = N×(N-1) ×(N-2) ×... ×2×1 通り
(m1,m2,...,mN)=(b(s1),b(s2),...,b(sN)) 割り付け方:Si=(s1,s2,...,sN) がN!通り

7 Si=(i1,i2,...,iN)という割り付け m1 = b(s1) m2 = b(s2) m1 = b(s2) ... ...
mN = b(sN)

8 Si=(i1,i2,...,iN)という割り付けを 観察する確率は？ P(m1=b(s1)) P(m2=b(s2)) P(m3=b(s3))
× P(m2=b(s2)) × P(m3=b(s3)) m1 = b(s1) m2 = b(s2) m1 = b(s2) × × × P(mN=b(sN)) ... ... mN = b(sN)

9 N!通りの確率 P(S1),P(S2),...,...,...,...,P(SN!) 最も大きなP(Si)となるSiは最尤推定割り付けがある

10 N!通りの確率 P(S1),P(S2),...,...,...,...,P(SN!) 計算する...? 1!=1 2!=2 3!=6 4!=24 5!=120 6!=720 7!=5040 8!=40320 9!=362880 10!=3,628,800 3百万 11!=39,916,800 12!=479,001,600 4.8億 15!=1.3 x 1012 20!=2.4 x 1018

11 多すぎてN!通りを計算できない N!通りを計算しないで、最尤割り付けがわかる？ 重みづけ最適化・重みづけマッチング
ハンガリアン・アルゴリズムなど

12 最尤推定割り付けがわかればそれが「答え」なのか？
P(Si) と P(Sj)とが第一位、第二位だとする P(Si)とP(Sj)とが等しかったら… P(Si)とP(Sj)とがほぼ等しかったら… P(Si)とP(Sj)とが数倍の違いしかなかったら…

13 第１,２,…n位割り付けがわかればそれが「答え」なのか？
「僅差」の割り付けがあったら、結局、どうしたらよいのかわからない 尤度が高い割り付けパターンを探すだけでは、解決しないかも

14 ある家族 「我が家の行方不明者mは、遺体b1,...,bNのうちのどれか１体だと言えますか？それとも、言いかねますか？」 ある遺体を保管しているところ 「この遺体bは、探されている行方不明者m1,...,mNのだれか１人だと言えますか？それとも、言いかねますか？」

15 N=3で考える ３ 体の遺体と３人の不明者を割り付ける場合の数 6通りの割り付け N=3 の場合 M1 M2 M3 仮説１ B1 B2 B3
仮説2 仮説3 仮説4 仮説5 仮説6 確率行列 B1 B2 B3 M1 P（M1=B1） P（M1=B2） M2 P（M2=B1） P（M2=B2） P（M2=B3） M3 P（M3=B1） P（M3=B2） P（M3=B3） 各仮説の尤度： ３つのペア全てでジェノタイプ が一致する確率

16 仮説2 6通りの割り付け N=3 の場合 M1 M2 M3 仮説１ B1 B2 B3 仮説2 仮説3 仮説4 仮説5 仮説6 確率行列 B1
P（M1=B1） P（M1=B2） P（M1=B3） M2 P（M2=B1） P（M2=B2） P（M2=B3） M3 P（M3=B1） P（M3=B2） P（M3=B3） この仮説の尤度＝P（M1=B1）xP（M2=B3）xP（M3=B2）

17 M1=B1,それ以外の割り付けはなんでもあり 仮説1+2
6通りの割り付け N=3　の場合 M1 M2 M3 仮説１ B1 B2 B3 仮説2 仮説3 仮説4 仮説5 仮説6 確率行列 B1 B2 B3 M1 P（M1=B1） P（M1=B2） P（M1=B3） M2 P（M2=B1） P（M2=B2） P（M2=B3） M3 P（M3=B1） P（M3=B2） P（M3=B3） この仮説の尤度＝P（M1=B1）xP（M2=B2）xP（M3=B3）

18 M1=B1,それ以外の割り付けはなんでもあり 仮説1+2
6通りの割り付け N=3　の場合 M1 M2 M3 仮説１ B1 B2 B3 仮説2 仮説3 仮説4 仮説5 仮説6 確率行列 B1 B2 B3 M1 P（M1=B1） P（M1=B2） P（M1=B3） M2 P（M2=B1） P（M2=B2） P（M2=B3） M3 P（M3=B1） P（M3=B2） P（M3=B3） この仮説の尤度＝P（M1=B1）xP（M2=B3）xP（M3=B2）

19 行列式(Determinant) 割り付けの場合ごとに掛け算をする「加える」要素と「引く」要素がある Wikipedia

20 パーマネント 割り付けの場合ごとに掛け算をする 全部を「加える」

21 パーマネント 行列式の計算は簡単で正確 この６通りの確率の和が３次正方行列のパーマネント パーマネントを求めるためのいくつかの方法がある
近似的 近似法を使うことでだいぶ速く計算できる。

22 パーマネントの計算方法 今日は、割愛 気になる方は
Wikipedia から情報の入手は可能です library(sets) library(gregmisc) # Ryser PermanentRyser<-function(A){ n<-length(A[,1]) p<-0 s<-as.set(1:n) for(i in 1:n){ scomb<-set_combn(s,i) #print(length(scomb)) c<-(-1)^i sumSameLen<-0 for(j in scomb){ #print(j) tmpprod<-1 for(k in 1:n){ tmpsum<-0 #print(k) for(l in j){ #print(A[k,l]) tmpsum<-tmpsum+A[k,l] } tmpprod<-tmpprod*tmpsum sumSameLen<-sumSameLen+tmpprod p<-p+c*sumSameLen p*(-1)^n # calculate based om difinition of permutation PermanentBeta<-function(A){ perms<-permutations(n,n) for(i in 1:length(perms[,1])){ t<-1 for(j in 1:n){ t<-t*A[j,perms[i,j]] p<-p+t p # Sinkhorn Balancing Scale sinkhorn<-function(A,epsilon= ){ B<-A scaleX<-scaleY<-rep(1,length(A[,1])) row_sum<-apply(A,1,sum) while(max(abs(t(row_sum)-1))>epsilon){ B<-B/row_sum scaleX<-scaleX/row_sum col_sum<-apply(B,2,sum) scaleY<-scaleY/col_sum B<-t(t(B)/col_sum) row_sum <- apply(B,1,sum) list(mat=B,scaleX=scaleX,scaleY=scaleY) # Huber/Law factor HuberLaw<-function(x){ x<-x x[which(x==0)]<-x[which(x==0)]+0.5 hl<-x mores<-which(hl>1) unMores<-which(hl<=1) hl[mores]<-hl[mores]+0.5*log(hl[mores])+exp(1)-1 hl[unMores]<-hl[unMores]*(exp(1)-1)+1 hl/(exp(1)) Minc<-function(x){ toosmall<-which(x<10^(-10)) x2<-x x2[toosmall]<-10^(-10) exp(lgamma(x2+1)*(1/x2)) # Approximation permanentEstimation<-function(A,iter=10000,epsilon=10^(-4),factor="hl"){ n<-length(A[1,]) shOut<-sinkhorn(A,epsilon) #print(shOut) B<-shOut[[1]] x<-shOut[[2]] y<-shOut[[3]] row_scale<-1/apply(B,1,max) C<-B*row_scale save_C<-C number_of_successes<-0 for(i in 1:iter){ column<-1 C<-save_C row_sums<-apply(C,1,sum) while(column<=n){ #print(column) if(factor=="hl"){ hl<-prod(HuberLaw(row_sums)) }else if(factor=="minc"){ hl<-prod(Minc(row_sums)) tmphl2<-HuberLaw(row_sums-C[,column]) tmphl2<-Minc(row_sums-C[,column]) hl2<-prod(tmphl2) #print(hl) #print(hl2) #print(tmphl2) row_probs<-C[,column]/tmphl2*hl2/hl #print(row_probs) # 重みづけサンプリング cumsumprob<-cumsum(row_probs) rrand<-runif(1) row_pick<-sum(cumsumprob<rrand)+1 #print(row_pick) if(row_pick==(n+1)){ column<-n+2 }else{ row_sums<-row_sums-C[,column] C[row_pick,]<-rep(0,n) column<-column+1 row_sums[row_pick]<-0 if(column==(n+1)){ number_of_successes=number_of_successes+1 hl_C<-prod(HuberLaw(row_sums)) hl_C<-prod(Minc(row_sums)) per_estimate<-hl_C*number_of_successes/iter per_estimate<-per_estimate/prod(row_scale)/prod(x)/prod(y) per_estimate ### MAIN FUNCTION ### # Calculate permanent CalcPermanent <- function(M,approx=T,epsilon=0.0001,iter=1000,factor="hl",beta=F, maxnfor.exact=16){ dimM <- dim(M) if(length(dimM)!=2) stop("M must have two dimensions") if(dimM[1]!=dimM[2]) stop("matrix must be a square matrix") n <- nrow(M) if(!approx){ if(n > maxnfor.exact) stop("n is too large for exact method. Use approx=T") else if(beta && n < 10) p <- PermanentBeta(M) else if(!beta || n >=10) p <- PermanentRyser }else if(approx){ if(any(M< 0)) stop("all elements of M must be non-negative values") p <- permanentEstimation(M,iter,epsilon,factor)

23 　近似法は速い 計算にかかる時間 （sec） 近似法 正確法 （N: 行列サイズ）

24 Fluctuation of estimation
もとがNｘN行列なら2ｘN 個の L(All) が出る 各列和・各行和 そのL(All) は推定値なので、誤差がある 最大値と最小値の差は0.102 注： scale 調整後

25 パーマネント計算 近似法の精度 「真のパーマネント」を大きな行列で計算するのは非現実的なので、「真の精度」を評価するのは難しいのですが
パーマネント計算　近似法の精度 「真のパーマネント」を大きな行列で計算するのは非現実的なので、「真の精度」を評価するのは難しいのですが Approximating the Permanent with Belief Propagation, by Bert Huang and Tony Jebara @ これは別なパーマネント近似法ですが…

26 尤度のNxN行列 「mi=bj、あとは何でもあり」に対応する(N-1)!仮説の確率を合算する この行列の各行の和は、どの行も等しい
各行の和は以下の和 「mi=b1、あとは何でもあり」 「mi=b2、あとは何でもあり」 … 「mi=bN、あとは何でもあり」 これは「miも何でもあり、他も何でもあり」だから

27 尤度のNxN行列 「mi=bj、あとは何でもあり」に対応する(N-1)!仮説の確率を合算する この行列の各列の和もやはり等しい
各列の和は、各行の和とも等しい

28 行を列に入れ替えても同じこと 家族が知りたいことにも、遺体保管者が知りたいことにも、答えられる
ある家族 「我が家の行方不明者mは、遺体b1,...,bNのうちのどれか１体だと言えますか？それとも、言いかねますか？」 列 ある遺体を保管しているところ 「この遺体bは、探されている行方不明者m1,...,mNのだれか１人だと言えますか？それとも、言いかねますか？」

29 尤度割合のNxN行列 尤度のNxN行列の各列の和、各行の和はすべて等しいので、その値で、尤度のNxN行列のすべての成分を割ってやる
各行、各列の和は、すべて１

30 ２つのNxN行列 確率行列 尤度割合行列 P(mi = bj) の行列
「m1=b2、あとは何でもあり」に対応する(N-1)!仮説の確率を合算

31 NxN の確率行列 割り付けの計算のために正方行列がほしい。 全ての遺体と行方不明者が１対１対応すると仮定
仮定できなければ、足りない分を一般集団で補う B1 B2 B3 ・・・・ BN M1 P（M1=B1） P（M1=B2） M2 P（M2=B1） P（M2=B2） P（M2=B3） M3 P（M3=B1） P（M3=B2） P（M3=B3） ： MN P（MN=B1） P（MN=B2） P（MN=B3）

32 NxN尤度割合行列 L(mi=bj) / L(ALL) の NxN 行列ができる 基準 v を満たしたペアの割り付けが決定する
・・・・ BN M1 0.9998 0.0000 M2 0.9953 M3 0.7 M4 0.0030 M5 0.0001 0.3 ： MN

33 行を列に入れ替えても同じこと 家族が知りたいことにも、遺体保管者が知りたいことにも、答えられる
ある家族 「我が家の行方不明者mは、遺体b1,...,bNのうちのどれか１体だと言えますか？それとも、言いかねますか？」 列 ある遺体を保管しているところ 「この遺体bは、探されている行方不明者m1,...,mNのだれか１人だと言えますか？それとも、言いかねますか？」

34 遺体引き取り　と　遺体引き渡し NxN尤度割合行列のセルの値を使って、「遺体引き取り」「遺体引き渡し」の判断ができるだろう 閾値は…

35 処理フロー Missings Bodies 血縁関係情報 遺体のＤＮＡ型 血縁者のＤＮＡ型 確率行列 パーマネント計算 尤度割合行列
引き取り・引き渡し判定

36 simulation data 想定した全体 手元のデータ 全ての行方不明者：104人 100家系 ＋ 重複4家系
100家系　＋　重複4家系 重複：行方不明者が複数いる家系 手元のデータ 100家系のうち68家系（incl. 重複３家系） 本当は全部使う予定だったが時間の関係で途中まで計算したところで割り付けをしたので家系も足りない状況になった。 104人の不明者のうち、54人をランダムに選択  家系も遺体も足りない状況

37 Simulation　結果

38 正解のペア全：30ペア v > 0.999 : 27 ペア ☆ v > 0.9 : 1 ペア ！ 家系として0.999 を満たす
全て正解 v > 0.9 : 1 ペア　！ 正解 家系として0.999 を満たす 1家系ー２遺体　が該当 0.999で感度　27/30　 ハズレなし

39 permanent計算後の確率（尤度比）が０でなかったペアの値をsort してplot したもの
正解ペア はずれペア 同一家系の2人 正解ペア

40 尤度比の自然対数 vs 1/(1-p) の自然対数 （p: 尤度割合）

41 可能性と課題 可能性 課題 「事前確率」を取り込める 「一般化」 申請のない行方不明者と発見されていない遺体のとりあつかい
性別・所持品等の情報を容易に取り込める 「一般化」 個人鑑定を同じ枠組みで考えることが(おそらく)可能 課題 申請のない行方不明者と発見されていない遺体のとりあつかい １家系に複数の行方不明者がいる場合 非独立な確率・・・行列を使った計算が苦手