■ - ryamadaの遺伝学・遺伝統計学メモ

複数のSNP型多型があり、それぞれ、0-1型形質のリスクを有するとする。
こちらの記事で示したように、各SNPのアレルには、アレル頻度に応じたリスクを持たせ、そのアレルに持たせたアレル特異的リスクの期待値は集団全体でゼロとなるように標準化するものとする。
個々のSNPのアレル特異的リスクを規定するパラメタについては、平均0、分散1の正規分布を仮定する。
このような条件で、HWE、diplotype的にもAdditive model,複数のローカスにおけるリスクもAdditive modelとしてやる。
SNP数と個人数とを引数に、このようなモデルの個人別形質リスクの分布とリスクの平均と分散とを求めてやるR関数を作ってみた。
こうして得られる個人別のリスクの分布において、発症閾値モデルでは、ある閾値以上の個人が発症する、というものである。
ローカス数、そのリスクを与える正規分布の条件などを検討するため。
ただし、閾値モデルにおいては、閾値の設定により、いずれにせよ、発症群の割合を決められるので、リスクを規定する正規分布条件については、あまり細かい条件はいらないか・・・

multiLocusRiskSum<-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/p1
u2<-(-w)/p2
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)
}

}

実行：" record<- multiLocusRiskSum(1000,1000,0.0001,0.995,"test4.txt")" 第１引数がSNP数、第２引数が個人人数、第３引数がsd、第４引数が、取り立てるquantile 閾値、第５引数が書き出しファイル名