■

フェノタイプの単位とジェノタイプの単位は違って、フェノタイプ=2とジェノタイプ=2の意味は違うのだけど、MSTで処理するときに、同じものとしてあつかわれているので、ちょっと微妙になっているかもしれない。困った。リスクの与え方を工夫しないといけない。フェノタイプが常に１以下がいいかなぁ・・・フェノタイプやリスクは相対的なものとしているので、値の範囲を決めても悪くはないはず。

１以下にするのは、ジェノタイプが0,1,2の値をとるので、関係がなさそうな感じに見える、というか、木のでき方が微妙にならない感じがするから。

リスクを＋のみでしらべるより、−のみで調べる方が大切だったかもしれない。

■ 

source("MST.R")
par(ask=F)

Np<-100
HN<-500
df<-2000
N<-280							#バリアントの数 i
ka<-1000
coa<-c(0,1000,2000,3000,4000,5000)
lcoa<-length(coa)
syu<-matrix(200,Np,lcoa)
syu2<-matrix(0,Np,lcoa)
colnames(syu)<-coa
nrvcv<-matrix(200,Np,lcoa)


afb<-0
for (i in 1:round(df/20)){
	afd<-(i-0.5)/df
	afb<-c(afb,rep(afd,round(N*(pexp(i/df,ka)-pexp((i-1)/df,ka)))))
}
raf<-afb[-1]
#caf<-c(0.02,0.02,0.03,0.04,0.05,0.06,0.08,0.10,0.20,0.30)
#af<-c(raf,caf)
af<-raf
hist(af,breaks=seq(0,max(af)+1/df,1/df),xlab="model allele frequency",main="allele frequency model")
N<-length(af)

for (inp in 1:Np){

ir<-al<-matrix(0,HN,N)		#各ローカスにおけるリスク、持っているアレル
P<-sr<-matrix(0,HN,lcoa)			#フェノタイプ、リスクの総和
iaf<-rep(0,N)		#真のアレル頻度
risk<-matrix(0,N,lcoa)

kijun<-sum(af)
plmi<-sign(rnorm(N, mean=0, sd=1))
for (icoa in 1:lcoa){
kr<-coa[icoa]
risk[,icoa]<-rnorm(N,mean=exp(-kr*af),sd=exp(-kr*af)*0.1)*plmi
kjn<-kijun/sum(abs(risk[,icoa])*af)
risk[,icoa]<-risk[,icoa]*kjn
plot(af,risk[,icoa],xlab="model allele frequency",main="allele frequency - risk")
}

#個人のアレルの有無振り分け iさんについて
for(i in 1:HN){
	a1<-a2<-tmp1<-tmp2<-rep(0,N)
	a1<-runif(N)					#ＲＶの有無を決定
	a2<-runif(N)
	tmp1[which(a1<af)]<-1
	tmp2[which(a2<af)]<-1
	al[i,]<-tmp1+tmp2				#アレルをいくつもつか
	for (icoa in 1:lcoa){
		sr[i,icoa]<-sum(risk[,icoa]*al[i,])		#リスク総和
	}
}


iaf<-apply(al,2,mean)/2			#アレル頻度確かめ
#plot(af,iaf,xlab="model allele frequency",ylab="real allele frequency")						#グラフを描画

for(icoa in 1:lcoa) plot(density(sr[,icoa]),xlab="sum(risk)",main="sum(risk) = Phenotype",xlim=c(-4,4),ylim=c(0,8))
P<-sr
syu2[inp,]<-apply(sr,2,var)
#sdfrac<-0.5			#遺伝率
#P<-sr*sdfrac+rnorm(N,sd=sqrt(var(sr))*(1-sdfrac))
#plot(sr,P)

#----------------ＭＳＴ・量的形質----------------


asdf<-apply(al,2,sum)
qwer<-which(asdf==0)
ala<-al[,-qwer]
jaf<-apply(ala,2,mean)/2

hist(jaf,breaks=seq(0,max(jaf)+1/HN,1/HN/2),xlab="real allele frequency",main="real allele frequency")

Nrv<-length(jaf[which(jaf<0.05)])
Ncv<-length(jaf[which(jaf>=0.05)])
for (icoa in 1:lcoa){
nrvcv[inp,icoa]<-Nrv
X<-cbind(P[,icoa],ala)
MST<-StatsMST(X,perm=TRUE,Nperm=100,tobeshuffled=c(1),dist.method="euclid")
plot(sort(MST$PermL),ylim=range(c(MST$L,MST$PermL)))
abline(h=MST$L,col=2)
title(paste(as.character(coa[icoa]),as.character(inp),as.character(Nrv),as.character(Ncv)))
ll<-MST$L
cat1<-seq(from=0,to=1,by=0.01)
cat2<-quantile(MST$PermL,cat1)
cat3<-cat2-ll
syu[inp,icoa]<-length(cat3[cat3<0])
}
}

boxplot(syu,xlab="kr (risk=rsign*exp(-kr*allelefrequency))")
syu

save.image("10051.RData")

なるべく、リスク総和の分布が似た形になると、検出力がリスクの与え方によって変わるのか、分布の仕方によって変わるのかの区別をつけられる、と考えた。

レアであるほどリスクが高くなり、その高くなる割合を変える、という前提は変えられないので、リスク総和が小さくなりがちなリスクの与え方になっているとき、リスクの値自体を増やすことにする。

sum(risk * allele frequency)の値が一緒になるようにすると、リスク総和の分布が近くなる。

■


1000人、N=200、観測されるRV150個くらい、krを変えたとき


500人、n=300、観測されるRV150個くらい、krを変えたとき

観測されるRVの数が揺れていて微妙なので、何人はこのRVを持ってる、と決めてしまって、誰が持っているかはランダムに決めるのがいいと思う。傾向を見たいだけなので、そこまで厳密であることが必要かわからないが。

■

十分な人数、十分なバリアントでリスクの与え方を変えたいのだが、人数を2倍にすると4倍くらいの時間がかかる気がする。どのくらいにするか難しいところ。

1000人、RV150個で、大体72秒。あるkrについて100回試行して2時間。

距離のはかり方について

距離のはかり方をなんとなくmanhattanからeuclidにしてみたら上手くいきそうな感じ。。。。

manhattanでCVをちょっと増やすとかより、こっちの違いについて考えるのが重要な感じがする。

euclidで色々な条件で試してみて、それで上手くいくならそれでいいし、上手くいかないときにmanhattanでやってみてうまくいったら、条件の違いを考えてみる。
→条件の違いで方法を変える。

プレゼンを作っているときに、もしかしたらEuclidのほうが説明しやすくて考えやすいかも、と思ったので、euclidでやってみよう。


横軸がkr

■

サンプルを200人とってきて、こういったRVセットについて考える、としたとき、どんなRVを選んで調べて見てもいいのだから、アレル頻度のところは実際は融通がききそうな気がしてきた。といっても、みつかるRVのアレル頻度は指数分布に従うだろうから、10-50%くらいのCVもあわせて少し探してもらうことにしてみる。

0-5%のRVのみ発生させたとき、リスクのあたえかたをかえると、

こんな感じ。何も分からない。

本当は、どんな頻度分布でも、RVにリスクがあれば検出されるようなのがいいのだけど、なぜか調査対象にRVしか存在しないと、上手くいかない。

論文のデータで、DHSについて、1%以下が120、5%以下が126、全体は129個というデータがあったから、その割合でやる。