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

• 2008/08/15にアレル特異的リスクの表現についてメモした(こちら)。
• 少し変える。変えると、２アレルのリスクの差がアレル頻度に依存しなくなる。
• 集団に0,1型の形質がある。そのprevalence をとする
• SNPがあり、0,1のアレルを持つ。それぞれのアレル頻度をとする。
• 集団での平均的リスクは、prevalence,で定まっており、SNPのアレル特異的リスクは、そこからのずれとして考えることとすると、２アレルのアレル特異的リスクをと置くことで、満足される。
  • がその条件である。
• このようなときの、集団の形質別・ジェノタイプ別の頻度は、HWEと、リスクのadditive modelを仮定すると以下のようになる。(に注意して式変形する)
  • 
  • 
  • さらに、フェノタイプ別のアレル本数の分布は
    • 
    • 
  • 書き換えると
    • 
    • 

• こう置いてもよい
  • 
  • 
  • さらに、フェノタイプ別のアレル本数の分布は
    • 
    • 
  • 書き換えると
    • 
    • 
• このアレル設定に対応させた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)
}