■
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) #RVの有無を決定 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) #----------------MST・量的形質---------------- 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)の値が一緒になるようにすると、リスク総和の分布が近くなる。