ジェノタイプで木を作ってフェノタイプを乗せる

#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)					#ＲＶの有無を決定
	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値の出し方が、今、パーミュテーションで出した統計量のうちもとの統計量より小さいものの比率なのだが、そしたらもっとパーミテーションしたほうがいいかも？