尤度とかそういうものについて その2

尤度について書いていたけどそれの続きのような違うような。
DNA鑑定についての講義があった。
担当教員に聞いてみればよかったのだが、質問したいことがまとまらないまま週末を迎えてしまったので、備忘録。
 
まあなんか事件があったと。
現場の遺留品からDNAを取って調べました。
現在は、マイクロサテライト(Single Tandem Repeat,STR)という、2〜5塩基の繰り返し配列を15ローカス調べている。
昔は9ローカスだったらしいので、この増量による精度の違いもまた後で考えることにして…
今、仮想的アレル頻度と検査法をここにおいておく。
・marker…調べるローカスの数。現在なら15か所。
・allele…1ローカスに存在するアレルの数。今は簡略化のため、各ローカスで1番から20番まであることにする。また、アレルの存在頻度はdirichlet関数で与える。

library(MCMCpack)
marker<- 15               # 調べるローカス数。各ローカスは独立と仮定する。
allele<- 20               # 存在するアレル数。
takes<- 10                # 現場で検出されるアレル数を規定する変数的な
A<- rdirichlet(marker,rep(1,allele))           # アレル存在頻度。行がローカス、列がアレル番号。
sus<- matrix(sample(1:allele,size=2*marker,replace=T),nc=2)     # 被疑者のもつアレルパターン。
                                                                # 父方と母方を区別していると思う。
B<- matrix(sample(1:allele,size=takes*marker,replace=T),nr=marker)  # 何をしているかというと…
C<- lapply((apply(B,1,unique)),sort)                # これが、現場で検出されたアレル。

ここで、結果が

 C
[[1]]  5  7  8 10 12 17 19 20
[[2]]  1  5  7 11 12 19
[[3]]  3  5  6  9 10 13 17 20
[[4]]  3  5  6 11 12 13 15 16
[[5]]  1  2  3 11 13 18 19 20
[[6]]  6  7  9 13 16 17 18 20
[[7]]  4  6  7  9 12 13 14
[[8]]  3  5  6  7  8  9 16 17 18 19
[[9]]  2  3  4  7  8  9 11 17 19
[[10]] 3  4  6  8  9 11 12 14 15
[[11]] 2  5  6  8 10 11 12 14 17 20
[[12]] 3  6  7  8  9 11 18 19
[[13]] 1  5  6  7 11 14 15 16
[[14]] 3  4  5  9 10 11 13 16
[[15]] 9 10 12 14 15 16
 sus
      [,1] [,2]
 [1,]    8   12
 [2,]    1   11
 [3,]    9   17
 [4,]    3    5
 [5,]   18    3
 [6,]    6   16
 [7,]    6    4
 [8,]    8    8
 [9,]    2    9
[10,]    4    3
[11,]   17    5
[12,]    6   19
[13,]    7   14
[14,]    9    4
[15,]   15   16

となったとする(ちょっといじった)。ここで、被疑者はすべてのローカスで、現場で検出されたアレルと同じアレルを持っている。
ここで、この被疑者の15ローカスのアレルの存在頻度は、
(父方染色体ローカス1のアレルの存在頻度)*(母方染色体ローカス1のアレルの存在頻度)

で、各ローカスは独立なので、各々かけ合わせる。
被疑者の15ローカスのアレルセットが存在する頻度を求める。

sus.prob<- 1
S<- list(cbind(1:marker,sus[,1]),cbind(1:marker,sus[,2]))
sus.prob<- prod(A[S[[1]]])*prod(A[S[[2]]])
sus.prob
[1] 1.471066e-45

となる。10^{-45}とかまあかぶることはありえないだろうから、データとしては
「この被疑者は現場にいたとして矛盾しない」
と出すらしいけど、そうしたら警察はウホッ逮捕となる。
 
が、ここからが質問事項で、逆に、上のCの結果を見ると、こんだけアレルあったら、まあ普通に誰でも適合すんじゃね?ということで、教員が自分のゲノムで試したら、なんと適合してしまったらしい。
ここで、現場で検出されたアレルから生まれうる15ローカスの組み合わせを考えると、各ローカスは独立なので

all.pattern<- (allele*allele)^marker                  # 全ローカス全アレルの組み合わせ総数。
variety<- 1
for(g in 1:length(C)){
  variety<- variety*(length(C[[g]])*length(C[[g]]))   # g番目のローカスについて、現場で検出されたアレルの組み合わせを考える。
  }
variety                                               # 現場から考えられるジェノタイプ総数。
[1] 1.172765e+27
variety/all.pattern                                   # 全体に対する存在確率。
[1] 1.092222e-12

と、兆単位になった。
これで、なんか数十人に一人かぶるからなんちゃらと言っていた気がするのだが…
 
参考資料
先生の法医学