本当のレアバリアントと見かけのレアバリアント

• 限定された本数の観測データがあって、標本におけるマイナーアレルの比率が低いとき、その中にはどれくらいの本当のレアバリアントが含まれ、どれくらいの「本当はコモンバリアントだけれど、たまたま少ない本数しかマイナーアレルを示さなかった」ものが含まれるかは、母集団のマイナーアレル頻度別バリアント分布に依存する話だが…
• そんな話→

# SNP,rare variantsの母集団でのアレル頻度別割合を
# 指数関数で近似する
# そのパラメタ
m <- 20
# アレル頻度別のバリアント割合をプロット
x <- seq(from=0,to=1,length=1000)
fx <- m * exp(-m*x)
plot(x,fx,type="l")

# N本の染色体をスクリーニングする
N <- 100

# 観察されるのは、N本のうち、0,1,2,...,N本が変異型として見つけられる
# 母集団の頻度に応じて２項分布になる
n <- 0:N
n2 <- N:0

# シミュレーションでやってみよう
# M箇所でそれをやれば
M <- 1000

ret <- rep(0,N+1)

tmp <- lgamma(N+1)-lgamma(n+1)-lgamma(n2+1)
for(i in 1:M){
	p <- rexp(1,m)
	tmp2 <- n*log(p)+n2*log(1-p)
	ret <- ret + exp(tmp+tmp2)
}
# すべてがワイルドタイプの場合n[1]を除いてプロットする
plot(n[-1],ret[-1])