■

#バリアントの数が60のときの、CNを変える

Date<-"09261"

source("MST.R")

Np<-100						#MST何回か
HN<-50000

df<-200
cnnu<-c(20,30,40,50,60,70)
Rnu<-length(cnnu)
syu<-matrix(0,Np,Rnu)	#boxplot用
colnames(syu)<-cnnu

for (nnp in 1:Rnu){
CN<-cnnu[nnp]
RN<-60
ka<-40
#アレル頻度
af<-afb<-0
for (i in 1:df){
	afd<-(i-1)/df
	if (afd==0) afd<-0.001
	afb<-c(afb,rep(afd,round(RN*(pexp(i/df,ka)-pexp((i-1)/df,ka)))))
}
af<-afb[-1]

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

caf<-crisk<-rep(0,CN)
caf<-runif(CN,0.2,0.5)

af<-c(af,caf)
risk<-c(risk,crisk)

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


eeee<-as.character(CN)
ssss<-paste(Date,"CN",eeee,".png")
png(ssss)          #デバイスドライバの用意
hist(af)
dev.off()

#個人のアレルの有無振り分け 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)                      #グラフを描画


#eeee<-as.character(df)
ssss<-paste(Date,"risk CN",eeee,".png")
png(ssss)          #デバイスドライバの用意
plot(density(sr))
dev.off()

#sdfrac<-0.98			#遺伝率
#P<-sr*sdfrac+rnorm(N,sd=sqrt(var(sr))*(1-sdfrac))
P<-sr    
#plot(sr,P)

#----------------ＭＳＴ・量的形質----------------
#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(CN*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])

}
save.image("09261vari.RData")
}

このafの与え方だと、合える頻度が低いものについては発生するが、高いものについてまったく発生しないので、少し考えないといけない。

多分CVをちょっといれるのと、afの与え方をよりそれらしくするのは、RNが40個や60個ぐらいのときは同じことなので、それでいこう。20個ぐらいのときはCVを入れる意味合いが変わってくる気がする。

pexpじゃなくてdexpでいい。。。dexpの逆関数があったらいいなぁ。。。