2009-08-22

構造化集団多次元空間

アレル頻度を多次元空間に正規分布を使って与える
個体は多次元空間の点においてＨＷＥを満足すると仮定する
サンプリングは多次元空間上において(正規)確率密度で行われる
異なるサンプリング集団間における集団構造の違いは、サンプリング確率密度分布のずれによって生じるものとする

#空間次元
K<-2
#マーカー数
Nx<-1000
#集団数
Np<-2

#マーカーのアレル頻度分布の平均値と分散
Mx<-matrix(runif(Nx*K),nrow=Nx)
Vx<-matrix(10*(runif(Nx*K))^2,nrow=Nx)
#マーカーのリスク係数
Rx<-rnorm(Nx)

#集団の分布の平均値と分散
#Mp<-matrix(runif(Np*K),nrow=Np)
Mp<-matrix(rnorm(Np*K,0,0.01),nrow=Np)
Vp<-matrix((runif(Np*K))^2,nrow=Np)
#Vp<-matrix(runif(Np*K),nrow=Np)

#集団別のサンプル数
Ns<-rep(0,Np)
Ns[1]<-100
Ns[2]<-100

#サンプルの位置を決める
Location<-function(ns,mp,vp){
ret<-matrix(rep(0,ns*length(mp)),nrow=ns)
for(i in 1:length(mp)){
ret[,i]<-rnorm(ns,mp[i],sqrt(vp[i]^2))
}
return(ret)
}

Ls1<-Location(Ns[1],Mp[1,],Vp[1,])
Ls2<-Location(Ns[2],Mp[2,],Vp[2,])

#サンプルの所在地でのマーカーのアレル頻度を決める
AlleleFreq<-function(L,M,V){
I<-length(L[,1])
J<-length(M[,1])
K<-length(L[1,])
ret<-matrix(rep(1,I*J),nrow=I)
for(i in 1:I){
for(j in 1:J){
for(k in 1:K){
 ret[i,j]<-ret[i,j]*pnorm(L[i,k],M[j,k],V[j,k])

}
}
}
return(ret)
}

AF1<-AlleleFreq(Ls1,Mx,Vx)
AF2<-AlleleFreq(Ls2,Mx,Vx)



#アレル頻度から、ディプロタイプをサンプリングする
SampleDiplotypeFromAlleleFreq<-function(af){
df1<-af*af
df2<-2*af*(1-af)
df3<-(1-af)*(1-af)
type<-c(0,1,2)
I<-length(af[,1])
J<-length(af[1,])
ret<-matrix(rep(0,I*J),nrow=I)
for(i in 1:I){
for(j in 1:J){
rrr<-runif(1)
if(rrr<df1[i,j]){
ret[i,j]<-0
}else if(rrr<df1[i,j]+df2[i,j]){
ret[i,j]<-1
}else{
ret[i,j]<-2
}
#ret[i,j]<-sample(type,1,replace=T,c(df1[i,j],df2[i,j],df3[1,j]))
}
}
return(ret)
}


GT1<-SampleDiplotypeFromAlleleFreq(AF1)
GT2<-SampleDiplotypeFromAlleleFreq(AF2)

#サンプルごとのリスクを計算する
R1<-GT1%*%Rx
R2<-GT2%*%Rx

#リスクのヒストグラムを描く
Hall<-hist(c(R1,R2))
hist(R1,breaks=Hall$breaks,ylim=c(0,max(Hall$counts)),col="blue",density=40,main="",xlab="",ylab="")
par(new=T)
hist(R2,breaks=Hall$breaks,ylim=c(0,max(Hall$counts)),col="red",density=40,main="",xlab="",ylab="")

#集団別のリスク平均
mean(R1)
mean(R2)

#リスクの累積分布を描く
maxrisk<-max(R1,R2)
minrisk<-min(R1,R2)
plot(sort(R1),ylim=c(minrisk,maxrisk),col="blue")
par(new=T)
plot(sort(R2),ylim=c(minrisk,maxrisk),col="red")

#サンプルの２次元位置分布を第一軸・第二軸で描く
max1<-max(Ls1[,1],Ls2[,1])
min1<-min(Ls1[,1],Ls2[,1])

max2<-max(Ls1[,2],Ls2[,2])
min2<-min(Ls1[,2],Ls2[,2])

plot(Ls1[,1],Ls1[,2],col="blue",xlim=c(min1,max1),ylim=c(min2,max2))
par(new =T)
plot(Ls2[,1],Ls2[,2],col="red",xlim=c(min1,max1),ylim=c(min2,max2))

#カウントの仕方(これは結局使わない)
count1<-apply(GT1,2,table)
count2<-apply(GT2,2,table)

#ケースコントロール集計

Rbind<-rbind(GT1,GT2)
phenotype<-c(rep(0,Ns[1]),rep(1,Ns[2]))

#それをトレンド検定
Trend<-function(Rbind,phenotype){
ret<-rep(0,length(Rbind[1,]))
for(i in 1:length(Rbind[1,])){
table<-table(Rbind[,i],phenotype)
w<-seq(from=0,to=length(table[,1])-1,by=1)
trendresult<-prop.trend.test(table[,1],table[,1]+table[,2],w)
ret[i]<-trendresult$p.value
}
return(ret)
}

trendP<-Trend(Rbind,phenotype)

#トレンド検定の累積分布を描く
plot(sort(trendP))