kr,v,pの3軸でplot
フェノタイプの分布を見たらkrが100以上だと分布がよくなかったので、krは0-100の範囲でやる。kr=100のとき、アレル頻度5%以上のCVのリスクは大体0になるので、kr=0:CDCV仮説、kr=100:MRV仮説、としてもいいだろう。
risk=k*rsign*exp(-kr*af)としたときに、kを10^(-9,-8.5,-8,,,,-6)と動かした時こんな感じで、知りたいところはこの図の左側のもやっとしているところと、その右の0に収束しているところとの間の部分あたりなので、もう少しkを刻んで分散が色々な値をとるようにする。今のkに、√1,√2,,,,,√9をかける感じにするとちょうどよさげ。kをr倍したら、分散は大体r^2になりそう。
そうすると、いま右の方で0に集まっているあたりはk=10^-6だから、もやもやっとしている左側は多分k=10^-6.5とか、もしかしたらk=10^-7。k=10^-8より小さいのとかは調べなくてもよかったかな・・・プログラムは回し始めてしまったので、明日の朝そこはちょっと変えよう。
ジェノタイプで木を作ってフェノタイプを乗せる
#tree-E,P-E,st source("MST.R") par(ask=F) Np<-100 #MST何回か HN<-200 df<-3500 N<-100 #バリアントの数 i ka<-600 kr<-c(0,10^(2:6)) lkr<-length(kr) syu<-matrix(1,Np,lkr) syul<-matrix(0,Np,lkr) colnames(syu)<-kr #allele frequency 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.03,0.03) af<-c(raf,caf) #hist(af,breaks=seq(0,max(af)+1/df,1/df),xlab="allele frequency",main="allele frequency model") N<-length(af) for (inp in 1:Np){ ir<-al<-matrix(0,HN,N) #各ローカスにおけるリスク、持っているアレル sr<-matrix(0,HN,lkr) wei<-rep(0,HN) P<-matrix(0,HN,lkr) #フェノタイプ、リスクの総和 iaf<-rep(0,N) #真のアレル頻度 risk<-matrix(0,N,lkr) plmi<-sign(rnorm(N, mean=0, sd=1)) for (i in 1:lkr) risk[,i]<-exp(-kr[i]*af)*plmi #for (i in 1:lkr) plot(af,risk[,i],xlab="model allele frequency",main="allele frequency - risk",xlim=c(0,0.01)) #個人のアレルの有無振り分け iさんについて for(i in 1:HN){ a1<-a2<-tmp1<-tmp2<-rep(0,N) a1<-runif(N) #RVの有無を決定 a2<-runif(N) tmp1[which(a1<af)]<-1 tmp2[which(a2<af)]<-1 al[i,]<-tmp1+tmp2 #アレルをいくつもつか for (j in 1:lkr) sr[i,j]<-sum(risk[,j]*al[i,]) #リスク総和 wei[i]<-sum(exp(-100*af)*al[i,]) } iaf<-apply(al,2,mean)/2 #アレル頻度確かめ #plot(af,iaf,xlab="model allele frequency",ylab="real allele frequency") #グラフを描画 #for (i in 1:lkr) plot(density(sr[,i]),xlab="sum(risk)",main="sum(risk) = Phenotype") #for (ikr in 1:lkr){ #sdfrac<-0.99 #遺伝率 #P[,ikr]<-sr[,ikr]*sdfrac+rnorm(HN,sd=sqrt(var(sr[,ikr]))*(1-sdfrac)) #P[,ikr]<-sr[,ikr] #P[,ikr]<-P[,ikr]/(abs(max(P[,ikr]))+abs(min(P[,ikr])))/2 #plot(sr[,ikr],P[,ikr,isc]) #} P<-sr StandardDist<-function(v,method="manhattan"){ N<-length(v) D<-dist(v,method=method) SumD<-sum(D) v/SumD*N*(N-1)/2 } asdf<-apply(al,2,sum) qwer<-which(asdf==0) ala<-al[,-qwer] alb<-apply(ala,2,StandardDist) jaf<-apply(ala,2,mean)/2 Nrv<-length(jaf[which(jaf<0.01)]) #hist(jaf,breaks=seq(0,max(jaf)+1/HN,1/HN/2),xlab="real allele frequency",main="real allele frequency") for (ikr in 1:lkr){ distM<-dist(alb,method="euclid") MST<-spantree(distM) di<-rep(0,HN-1) for (i in 1:HN-1){ di[i]<-sqrt((P[i+1,ikr]-P[MST[[1]][i],ikr])^2+(MST[[2]][i])^2) } L<-sum(di) PermL<-rep(0,100) for (nper in 1:100){ t<-sample(1:HN) PermP<-P[,ikr] for (i in 1:HN) PermP[i]<-P[t[i],ikr] di<-rep(0,HN-1) for (i in 1:HN-1){ di[i]<-sqrt((PermP[i+1]-PermP[MST[[1]][i]])^2+MST[[2]][i]^2) } PermL[nper]<-sum(di) } plot(sort(PermL),ylim=range(c(L,PermL))) abline(h=L,col=2) title(paste(as.character(inp),as.character(kr[ikr]),as.character(Nrv))) ll<-L cat1<-seq(from=0,to=1,by=0.01) cat2<-quantile(PermL,cat1) cat3<-cat2-ll syu[inp,ikr]<-length(cat3[cat3<0])/100 syul[inp,ikr]<-L } save.image("10132.RData") } colnames(syu)<-kr boxplot(syu,xlab="kr (risk=rsign*exp(-kr*allelefrequency))",ylab="p値",ylim=c(0,0.8)) apply(syu,2,var) title("Euclid,st")
euclid距離、標準化ありのとき
manhattan距離のときは
di[i]<-abs(P[i+1,ikr]-P[MST1[i],ikr])になるくらい。結果は計算中だけど、標準化したほうがよさげ。
p値の出し方が、今、パーミュテーションで出した統計量のうちもとの統計量より小さいものの比率なのだが、そしたらもっとパーミテーションしたほうがいいかも?
パラメタ
人数は十分な人数で固定して、多型の箇所数とアレル頻度は実際のデータに似せたものに固定する。
パーミテーションの回数は分布を出したいだけなので固定
フェノタイプの幅については、パラメタというり、検定手法の一部みたいなところがあるので、これも固定する。
リスクの与え方と遺伝率が問題。
・遺伝率が高いときにどういったリスクの与え方でワークするか
・ワークするリスクの与え方について、遺伝率を変えるとどうなるか
を調べるのがよいと思われる。
でも関連性あるかもしれないから同時に調べるべきか。。。
評価
パーミテーションテストでp値を出して、p<0.01であるものが99%以上のとき、検定が上手くいっていることにする。MSTの辺の長さの和が短い方がいいから、片側検定になると思う。まずp値の出し方を調べる。
p<0.01のものが何%か、というデータとして出す方が適当な気がしてきた。そうしよう。100%ならとてもいいわけだし。
3dplot
al<-matrix(0,36,2) al[,1]<-c(rep(2,9),rep(1,12),rep(0,15)) al[,2]<-c(2,2,1,1,1,0,0,0,0,2,2,2,1,1,1,1,0,0,0,0,0,2,2,2,2,1,1,1,1,1,0,0,0,0,0,0) sr<-al[,1]*1+al[,2]*1 P<-sr+rnorm(36,0,0.1) rgl.clear() rgl.bg(col="white") x<-cbind(P,al) d<-dist(x) stree<-spantree(d) library(rgl) plot3d(x) t<-matrix(0,length(x[,1]),2) for (i in 2:length(x[,1])){ ki<-c(i,as.numeric(stree$kid[i-1])) #t[,i]<-c(i,stree$kid[i]) rgl.lines(x[ki,1],x[ki,2],x[ki,3],col="red") } tmp1<-sample(1:36) tmp2<-sample(1:36) asl<-al for(i in 1:36) al[i,1]<-asl[tmp1[i],1] for(i in 1:36) al[i,2]<-asl[tmp2[i],2] rgl.clear() rgl.bg(col="white") x<-cbind(P,al) d<-dist(x) stree<-spantree(d) library(rgl) plot3d(x) t<-matrix(0,length(x[,1]),2) for (i in 2:length(x[,1])){ ki<-c(i,as.numeric(stree$kid[i-1])) #t[,i]<-c(i,stree$kid[i]) rgl.lines(x[ki,1],x[ki,2],x[ki,3],col="red") }
最初に関連度高いときのMST描いて、次にそれをパーミテーションしたのを描く par(ask=T)がきかなかった。
意外と3次元だと長さ分かりにくい。