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

多人数一括ＤＮＡプロファイリングのための
確率計算法に関する考察
法数学勉強会
2011/09/10
京都大学大学院医学研究科
統計遺伝学分野
奈良原舞子山田亮
状況
• 大災害が発生
• 多数の行方不明者
• 多数の身元不明遺体
• 外見や所持品などから身元が特定された遺
体はすでに遺族に返還されている。
• 残っている遺体の手がかりは、主に遺伝情報
使えるデータ
• 遺体
– ジェノタイプ
– 発見時の状況など
• 行方不明者
– 家系情報
– 家族のジェノタイプ
個人の鑑定
簡単に描くと
行方不明者
Missing
m1
身元不明者
found Body
b1
集団の誰か
集団の鑑定
行方不明者
Missing
m1
b1
m2
b2
m3
b3
..
..
..
..
..
..
mN
bN
身元不明者
found Body
集団の鑑定
N人をN体に割り付ける
• N! = N×(N-1) ×(N-2) ×... ×2×1 通り
– (m1,m2,...,mN)=(b(s1),b(s2),...,b(sN))
– 割り付け方:Si=(s1,s2,...,sN) がN!通り
Si=(i1,i2,...,iN)という割り付け
m1 = b(s1)
...
m2 = b(s2)
m1 = b(s2)
...
mN = b(sN)
Si=(i1,i2,...,iN)という割り付けを
観察する確率は？
P(m1=b(s1))× P(m2=b(s2))× P(m3=b(s3))
m1 = b(s1)
...
m2 = b(s2)
×
...
m1 = b(s2)
×
× P(mN=b(sN))
mN = b(sN)
N!通りの確率
• P(S1),P(S2),...,...,...,...,P(SN!)
• 最も大きなP(Si)となるSiは最尤推定割り付け
がある
• N!通りの割り付けのすべての確率を計算して、
最大の場合を見つける？
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
多すぎてN!通りを計算できない
• N!通りを計算しないで、最尤割り付けがわか
る？
– 重みづけ最適化・重みづけマッチング
• ハンガリアン・アルゴリズムなど
最尤推定割り付けがわかればそれが
「答え」なのか？
• P(Si) と P(Sj)とが第一位、第二位だとする
– P(Si)とP(Sj)とが等しかったら…
– P(Si)とP(Sj)とがほぼ等しかったら…
– P(Si)とP(Sj)とが数倍の違いしかなかったら…
第１,２,…n位割り付けがわかればそれ
が「答え」なのか？
• 「僅差」の割り付けがあったら、結局、どうした
らよいのかわからない
尤度が高い割り付けパター
ンを探すだけでは、解決しな
いかも
• ある家族
– 「我が家の行方不明者mは、遺体b1,...,bNのうち
のどれか１体だと言えますか？それとも、言いか
ねますか？」
• ある遺体を保管しているところ
– 「この遺体bは、探されている行方不明者
m1,...,mNのだれか１人だと言えますか？それと
も、言いかねますか？」
N=3で考える
• ３体の遺体と３人の不明者を割り付ける場合
の数
6通りの割り付け
N=3 の場合
M1
M2
M3
仮説１
B1
B2
B3
確率行列
B1
B2
B3
M1
P（M1=B1）
P（M1=B2）
P（M1=B2）
仮説2
B1
B3
B2
M2
P（M2=B1）
P（M2=B2）
P（M2=B3）
仮説3
B2
B1
B3
M3
P（M3=B1）
P（M3=B2）
P（M3=B3）
仮説4
B2
B3
B1
仮説5
B3
B1
B2
仮説6
B3
B2
B1
各仮説の尤度：
３つのペア全てでジェノタイプが一致する確率
仮説2
6通りの割り付け
N=3 の場合
M1
M2
M3
仮説１
B1
B2
B3
確率行列
B1
B2
B3
M1
P（M1=B1）
P（M1=B2）
P（M1=B3）
仮説2
B1
B3
B2
M2
P（M2=B1）
P（M2=B2）
P（M2=B3）
仮説3
B2
B1
B3
M3
P（M3=B1）
P（M3=B2）
P（M3=B3）
仮説4
B2
B3
B1
仮説5
B3
B1
B2
仮説6
B3
B2
B1
この仮説の尤度＝P（M1=B1）xP（M2=B3）xP（M3=B2）
M1=B1,それ以外の割り付けはなんでもあり
仮説1+2
6通りの割り付け
N=3 の場合
M1
M2
M3
仮説１
B1
B2
B3
確率行列
B1
B2
B3
M1
P（M1=B1）
P（M1=B2）
P（M1=B3）
仮説2
B1
B3
B2
M2
P（M2=B1）
P（M2=B2）
P（M2=B3）
仮説3
B2
B1
B3
M3
P（M3=B1）
P（M3=B2）
P（M3=B3）
仮説4
B2
B3
B1
仮説5
B3
B1
B2
仮説6
B3
B2
B1
この仮説の尤度＝P（M1=B1）xP（M2=B2）xP（M3=B3）
M1=B1,それ以外の割り付けはなんでもあり
仮説1+2
6通りの割り付け
N=3 の場合
M1
M2
M3
仮説１
B1
B2
B3
確率行列
B1
B2
B3
M1
P（M1=B1）
P（M1=B2）
P（M1=B3）
仮説2
B1
B3
B2
M2
P（M2=B1）
P（M2=B2）
P（M2=B3）
仮説3
B2
B1
B3
M3
P（M3=B1）
P（M3=B2）
P（M3=B3）
仮説4
B2
B3
B1
仮説5
B3
B1
B2
仮説6
B3
B2
B1
この仮説の尤度＝P（M1=B1）xP（M2=B3）xP（M3=B2）
行列式(Determinant)
割り付けの場合ごとに掛け算をする
「加える」要素と「引く」要素がある
Wikipedia
パーマネント
割り付けの場合ごとに掛け算をする
全部を「加える」
パーマネント
• 行列式の計算は簡単で正確
• この６通りの確率の和が３次正方行列のパー
マネント
• パーマネントを求めるためのいくつかの方法
がある
– 正確
– 近似的
• 近似法を使うことでだいぶ速く計算できる。
パーマネントの計算方法
• 今日は、割愛
– 気になる方は
• Wikipedia http://en.wikipedia.org/wiki/Permanent か
ら情報の入手は可能です
近似法は速い
（sec）
計算にかかる時間
正確法
近似法
（N: 行列サイズ）
Fluctuation of estimation
• もとがNｘN行列なら2ｘN 個の L(All) が出る
– 各列和・各行和
• そのL(All) は推定値なので、誤差がある
– 最大値と最小値の差は0.102
注：
scale 調整後
パーマネント計算近似法の精度
• 「真のパーマネント」を大きな行列で計算する
のは非現実的なので、「真の精度」を評価す
るのは難しいのですが
Approximating the Permanent with Belief
Propagation, by Bert Huang and Tony Jebara
@ http://www.cs.columbia.edu/~bert/permanentTR.pdf
これは別なパーマネント近似法ですが…
尤度のNxN行列
• 「mi=bj、あとは何でもあり」に対応する(N-1)!
仮説の確率を合算する
• この行列の各行の和は、どの行も等しい
– 各行の和は以下の和
•
•
•
•
「mi=b1、あとは何でもあり」
「mi=b2、あとは何でもあり」
…
「mi=bN、あとは何でもあり」
– これは「miも何でもあり、他も何でもあり」だから
尤度のNxN行列
• 「mi=bj、あとは何でもあり」に対応する(N-1)!
仮説の確率を合算する
• この行列の各列の和もやはり等しい
• 各列の和は、各行の和とも等しい
行を列に入れ替えても同じこと
家族が知りたいことにも、遺体保管者
が知りたいことにも、答えられる
• 行
– ある家族
• 「我が家の行方不明者mは、遺体b1,...,bNのうちのどれか１
体だと言えますか？それとも、言いかねますか？」
• 列
– ある遺体を保管しているところ
• 「この遺体bは、探されている行方不明者m1,...,mNのだれ
か１人だと言えますか？それとも、言いかねますか？」
尤度割合のNxN行列
• 尤度のNxN行列の各列の和、各行の和はす
べて等しいので、その値で、尤度のNxN行列
のすべての成分を割ってやる
• それを「尤度割合のNxN行列」とする
– 各行、各列の和は、すべて１
２つのNxN行列
• 確率行列
– P(mi = bj) の行列
• 尤度割合行列
– 「m1=b2、あとは何でもあり」に対応する(N-1)!仮
説の確率を合算
NxN の確率行列
• 割り付けの計算のために正方行列がほしい。
• 全ての遺体と行方不明者が１対１対応すると仮定
• 仮定できなければ、足りない分を一般集団で補う
B1
B2
B3
M1
P（M1=B1）
P（M1=B2）
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）
・・・・
・・・・
・・・・
・・・・
・・・・
BN
NxN尤度割合行列
• L(mi=bj) / L(ALL) の NxN 行列ができる
• 基準 v を満たしたペアの割り付けが決定する
B1
B2
B3
– v0.9998
= 0.99
としたら、
0.0000 0.0000
M1
M2
0.0000 0.9953 0.0000
M3
0.0000 0.0000 0.7
M4
0.0000 0.0030 0.0000
M5
0.0001 0.0000 0.3
：
：
：
：
：
：
：
：
：
：
：
：
MN
0.0000 0.0000 0.0000
B4
B5
・・・・
・・・・
・・・・
BN
行を列に入れ替えても同じこと
家族が知りたいことにも、遺体保管者
が知りたいことにも、答えられる
• 行
– ある家族
• 「我が家の行方不明者mは、遺体b1,...,bNのうちのどれか１
体だと言えますか？それとも、言いかねますか？」
• 列
– ある遺体を保管しているところ
• 「この遺体bは、探されている行方不明者m1,...,mNのだれ
か１人だと言えますか？それとも、言いかねますか？」
遺体引き取りと遺体引き渡し
• NxN尤度割合行列のセルの値を使って、「遺
体引き取り」「遺体引き渡し」の判断ができる
だろう
• 閾値は…
処理フロー
Missings
血縁関係情報
血縁者のＤＮＡ型
Bodies
遺体のＤＮＡ型
確率行列
パーマネント計算
尤度割合行列
引き取り・引き渡し判定
simulation data
• 想定した全体
– 全ての行方不明者：104人
– 100家系＋重複4家系
• 重複：行方不明者が複数いる家系
• 手元のデータ
– 100家系のうち68家系（incl. 重複３家系）
• 本当は全部使う予定だったが時間の関係で途中まで
計算したところで割り付けをしたので家系も足りない状
況になった。
– 104人の不明者のうち、54人をランダムに選択
–  家系も遺体も足りない状況
Simulation 結果
• 正解のペア全：30ペア
– v > 0.999 : 27 ペア ☆
• 全て正解
– v > 0.9 : 1 ペア！
• 正解
– 家系として0.999 を満たす
• 1家系ー２遺体が該当
– 正解
0.999で感度 27/30
ハズレなし
• permanent計算後の確率（尤度比）が０でな
かったペアの値をsort してplot したもの
正解ペア
はずれペア
同一家系の2人
正解ペア
尤度比の自然対数 vs 1/(1-p) の自然対数（p: 尤度割合）
可能性と課題
• 可能性
– 「事前確率」を取り込める
• 性別・所持品等の情報を容易に取り込める
– 「一般化」
• 個人鑑定を同じ枠組みで考えることが(おそらく)可能
• 課題
– 申請のない行方不明者と発見されていない遺体のと
りあつかい
– １家系に複数の行方不明者がいる場合
• 非独立な確率・・・行列を使った計算が苦手

Download Report