アレル特異的リスクの表現その2

  • 2008/08/15にアレル特異的リスクの表現についてメモした(こちら)。
  • 少し変える。変えると、2アレルのリスクの差がアレル頻度に依存しなくなる。
  • 集団に0,1型の形質がある。そのprevalence をqとする
  • SNPがあり、0,1のアレルを持つ。それぞれのアレル頻度をp_0,p_1とする。
  • 集団での平均的リスクは、prevalence,qで定まっており、SNPのアレル特異的リスクは、そこからのずれとして考えることとすると、2アレルのアレル特異的リスクu_0,u_1u_0=p_1w,u_1=-p_0wと置くことで、満足される。
    • p_0u_0+p_1u_1=0がその条件である。
  • このようなときの、集団の形質別・ジェノタイプ別の頻度は、HWEと、リスクのadditive modelを仮定すると以下のようになる。(p_0u_0+p_1u_1=0に注意して式変形する)
    • p_0^2q(1+2u_0),2p_0p_1q(1+u_0+u_1),p_1^2q(1+2u_1)
    • p_0^2(1-q(1+2u_0)),2p_0p_1(1-q(1+u_0+u_1)),p_1^2(1-q(1+2u_1))
    • さらに、フェノタイプ別のアレル本数の分布は
      • 2p_0q(1+u_0),2p_1q(1+u_1)
      • 2p_0(1-q(1+u_0)),2p_1(1-q(1+u_1))
    • 書き換えると
      • 2p_0q+2p_0p_1qw,2p_1q-2p_0p_1qw
      • 2p_0(1-q)-2p_0p_1qw,2p_1(1-q)+2p_0p_1qw
  • こう置いてもよい
    • p_0^2(q+2u_0),2p_0p_1(q+u_0+u_1),p_1^2(q+2u_1)
    • p_0^2(1-(q+2u_0)),2p_0p_1(1-(q+u_0+u_1)),p_1^2(1-(q+2u_1))
    • さらに、フェノタイプ別のアレル本数の分布は
      • 2p_0(q+u_0),2p_1(q+u_1)
      • 2p_0(1-(q+u_0)),2p_1(1-(q+u_1))
    • 書き換えると
      • 2p_0q+2p_0p_1w,2p_1q-2p_0p_1w
      • 2p_0(1-q)-2p_0p_1w,2p_1(1-q)+2p_0p_1w
  • このアレル設定に対応させたRの関数は以下の通り
multiLocusRiskSum2<-function(nummarker,numsample,s,quant,out="out.txt"){
k<-nummarker
numsample<-numsample
sd<-s
p1<-runif(n=k)
p2<-1-p1
w<-rnorm(n=k,sd=sd)
u1<-w*p2
u2<-(-w)*p1
x<-matrix(runif(n=k*numsample),nrow=k)
x2<-matrix(runif(n=k*numsample),nrow=k)
bool1<-(x<=p1)
bool2<-(x>p1)
bool3<-(x2<=p1)
bool4<-(x2>p1)
z<-t(bool1)%*%u1+t(bool2)%*%u2+t(bool3)%*%u1+t(bool4)%*%u2
#z<-rowSums(y)
plot(density(z))
ret<-list(NULL)
ret_no<-0
ret[[ret_no<-ret_no+1]]<-list(c(c("length","mean","var"),c(length(z),mean(z),var(z))))
quantile<-quantile(z,quant)
ret[[ret_no<-ret_no+1]]<-c("quantile",quantile)
data<-c(c(0,0),c(p1,p1))
data<-c(data,c(0,0),c(u1,u1))
data<-c(data,c(0,0),c(u2,u2))
counter<-0
for(i in 1:length(z)){
 if(z[i]>quantile){
  data<-c(data,c(i,z[i]),c(sign((t(bool1)[i,]-0.5)*u1),sign((t(bool3)[i,]-0.5)*u1)))
  #data<-c(data,c(i,z[i]),c(t(bool1)[i,],t(bool3)[i,]))
  counter<-counter+1
  #ret[[ret_no<-ret_no+1]]<-c(c(i,z[i]),c(bool1[i,],bool3[i,]))
 }
}
numcol=2+k*2
matrixdata<-matrix(data,ncol=numcol,byrow=TRUE)
ret[[ret_no<-ret_no+1]]<-t(matrixdata)
write.table(ret[3],file=out,sep="\t")
return(ret)
}

multiLocusRiskSum2<-function(nummarker,numsample,s,quant,out="out.txt"){
k<-nummarker
numsample<-numsample
sd<-s
p1<-runif(n=k)
p2<-1-p1
w<-rnorm(n=k,sd=sd)
u1<-w*p2
u2<-(-w)*p1
x<-matrix(runif(n=k*numsample),nrow=k)
x2<-matrix(runif(n=k*numsample),nrow=k)
bool1<-(x<=p1)
bool2<-(x>p1)
bool3<-(x2<=p1)
bool4<-(x2>p1)
z<-t(bool1)%*%u1+t(bool2)%*%u2+t(bool3)%*%u1+t(bool4)%*%u2
#z<-rowSums(y)
plot(density(z))
ret<-list(NULL)
ret_no<-0
ret[[ret_no<-ret_no+1]]<-list(c(c("length","mean","var"),c(length(z),mean(z),var(z))))
quantile<-quantile(z,quant)
ret[[ret_no<-ret_no+1]]<-c("quantile",quantile)
data<-c(c(0,0),c(p1,p1))
data<-c(data,c(0,0),c(w,w))
data<-c(data,c(0,0),c(u1,u1))
data<-c(data,c(0,0),c(u2,u2))
counter<-0
value<-p1*p2*w*p1*p2*w
data<-c(data,c(0,0),c(value,value))
for(i in 1:length(z)){
 if(z[i]>quantile){
  data<-c(data,c(i,z[i]),c(sign((t(bool1)[i,]-0.5)*u1),sign((t(bool3)[i,]-0.5)*u1)))
  #data<-c(data,c(i,z[i]),c(t(bool1)[i,],t(bool3)[i,]))
  counter<-counter+1
  #ret[[ret_no<-ret_no+1]]<-c(c(i,z[i]),c(bool1[i,],bool3[i,]))
 }
}
numcol=2+k*2
matrixdata<-matrix(data,ncol=numcol,byrow=TRUE)
matrixdata<-t(matrixdata)
matrixdata<-matrixdata[,order(matrixdata[2,])]
matrixdata<-matrixdata[order(matrixdata[,2]),]
ret[[ret_no<-ret_no+1]]<-matrixdata


write.table(ret[3],file=out,sep="\t")
return(ret)
}