バリアントの数を変えたとき

source("MST.R")
Np<-100							#MST何回か
HN<-10000
syu<-matrix(0,Np,7)	#boxplot用

nnu<-c(20,40,60,80,100,120,160)

for (nnp in 1:7){
N<-nnu[nnp]

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

if (nnp==1) af<-c(rep(0.001,7),rep(0.01,5),rep(0.02,3),rep(0.03,2),rep(0.04,1),rep(0.05,1),0.07)
if (nnp==2) af<-c(rep(0.001,14),rep(0.01,9),rep(0.02,6),rep(0.03,3),rep(0.04,3),rep(0.05,2),0.06,0.07,0.10)
if (nnp==3) af<-c(rep(0.001,20),rep(0.01,14),rep(0.02,10),rep(0.03,7),rep(0.04,4),rep(0.05,2),0.06,0.07,0.10)
if (nnp==4) af<-c(rep(0.001,28),rep(0.01,18),rep(0.02,12),rep(0.03,8),rep(0.04,5),rep(0.05,4),0.06,0.06,0.07,0.07,0.10)
if (nnp==5) af<-c(rep(0.001,35),rep(0.01,24),rep(0.02,16),rep(0.03,10),rep(0.04,7),rep(0.05,4),0.06,0.06,0.07,0.10)
if (nnp==6) af<-c(rep(0.001,41),rep(0.01,28),rep(0.02,19),rep(0.03,13),rep(0.04,8),rep(0.05,5),0.06,0.06,0.06,0.07,0.07,0.10)
if (nnp==7) af<-c(rep(0.001,55),rep(0.01,37),rep(0.02,24),rep(0.03,16),rep(0.04,11),rep(0.05,8),0.06,0.06,0.06,0.06,0.07,0.07,0.07,0.10,0.10)

kr<-40
risk<-exp(-kr*af)*sign(rnorm(N, mean=0, sd=1))

#個人のアレルの有無振り分け 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				#アレルをいくつもつか
	sr[i]<-sum(risk*al[i,])		#リスク総和
}

#for (i in 1:N) iaf[i]<-mean(al[,i])/2	#アレル頻度確かめ
#plot(af,iaf)                      #グラフを描画
#sdfrac<-0.98			#遺伝率
#P<-sr*sdfrac+rnorm(N,sd=sqrt(var(sr))*(1-sdfrac))
P<-sr    


#----------------ＭＳＴ・量的形質----------------
#par(ask=FALSE)

Ncaco<-200

for (counte in 1:Np){
	Cmst<-sample(1:HN,Ncaco,replace=FALSE)
	X<-cbind(P[c(Cmst)],al[c(Cmst),])
	
	MST<-StatsMST(X,perm=TRUE,Nperm=100,tobeshuffled=c(1),dist.method="manhattan")
	plot(sort(MST$PermL),ylim=range(c(MST$L,MST$PermL)))
	abline(h=MST$L,col=2)
	title(N*1000+counte)
	ll<-MST$L
	cat1<-seq(from=0,to=1,by=0.01)
	cat2<-quantile(MST$PermL,cat1)
	cat3<-cat2-ll
	syu[counte,nnp]<-length(cat3[cat3<0])

}
}
colnames(syu)<-nnu