検出したいものがマスクされる

• 見たいアレルがプールされたサンプル中にあるかどうかを調べたい
• 実験のコールミスがある
• 標本中の真のアレル別本数がわかっているものとする
• そこから指定本数を取り出して実験する
• 取り出された分の真のアレル本数分布が出る
• それに、さらにコールエラーを加味して計算したものが、観測されるアレル別本数の分布

# N本の染色体
N<-100
# mutant alleleの本数
a<-3
# wild alleleの本数
A<-N-a

# m本を取り出したら、そのうち、i本がmutant alleleである確率

prob.i<-function(N,a,m,i){
	four<-c(i,a-i,m-i,N-(a+m)+i)
	ret<-0
	if(prod(four>=0)){
		ret<-exp(lgamma(N-a+1)+lgamma(a+1)+lgamma(m+1)+lgamma(N-m+1)-lgamma(N+1)-sum(lgamma(four+1)))
	}
	ret
}
m<-3
prob.i(N,a,m,1)

# アレルAであるときに、観測アレルがA,a,non-(A/a)である確率がPであるとしたとき
# k本のアレルAを観測して、x1,x2,x3なる本数で観測する確率は
# dmultinom(x,prob=P)

# アレルa/Aがk1,k2本であるときに、P,Qなる観測確率であるとき、
# x1,x2,x3なる本数で観測する確率
dmultinom.SNP<-function(x,a,prob1=P,prob2=Q){
	N<-sum(a)
	if(sum(x)!=N){
		return(0)
	}
	ret<-0
	for(i in 0:min(a[1],x[1])){
		for(j in 0:min(a[1]-i,x[2])){
			six<-c(i,j,a[1]-(i+j),0,0,0)
			six[4:6]<-x-six[1:3]
			if(prod(six>=0)){
				tmp<-dmultinom(six[1:3],prob=prob1,log=TRUE)+dmultinom(six[4:6],prob=prob2,log=TRUE)
				ret<-ret+exp(tmp)
			}
		}
	}
	ret
}

dmultinom.SNP(c(1,1000,0),c(1,1000),prob1=c(0.99,0.01,0),prob2=c(0.01,0.99,0.0))
#########
N<-100
a<-1
m<-50
k<-0:1
out.prob.i<-rep(0,length(k))
for(i in 1:length(k)){
	out.prob.i[i]<-prob.i(N,a,m,k[i])
}

P1<-c(0.99,0.01/3,0.01*2/3)
P2<-c(0.01/3,0.99,0.01*2/3)
ret2<-c()
for(i in 1:length(k)){
	A<-c(k[i],m-k[i])
	ret2<-c(ret2,out.prob.i[i]*dmultinom.SNP(c(a,m-a,0),A,prob1=P1,prob2=P2))
}
ret2