2009-07-24

■

先日、複合遺伝性疾患のリスクをヒストグラム表示するための簡易Rコマンドを書いた。
その続き。
今、ある複合遺伝性疾患があって、その発病リスクは環境要因と遺伝要因とで構成され、それぞれは相互に独立であり、また、それぞれは正規分布を仮定できるとする。
今、遺伝要因を構成する相互に独立な多座位のうち、ある１つの座位のリスクに着目する。
以下に示すRの関数の説明は次の通り。
シミュレーション回数N　デフォルトは10000
疾患の有病率prevalence
有病率のうち、この試行のモデルで説明できる分が1-phenocopy, デフォルトはphenocopy=0
遺伝率（環境要因と遺伝要因を足した、全分散のうち、遺伝要因が占める割合）heritability
遺伝要因のうち、着目している座位が担う割合 locusFrac(1ならこの座位のみということ）（厳密には、着目座位以外のリスクに正規分布を仮定し、着目座位のそれは、２アレルの作る３ディプロタイプでのそれになるので、異なってくるが。
着目座位のリスクアレルの頻度alleleFreq
３ディプロタイプの頻度がHWEからずれていてもよいようにその因子F デフォルトは0でそのときはHWE
homo1, hetero, homo2のリスクは線形をデフォルト(k=0)とするが、k=1にすれは、heteroのリスクはhomo1と同じになり、k=-1とすればheteroのリスクはhomo2と同じになる。そのためのパラメタがk

掲載図は、以下の実行例の最初の方のヒストグラム（３色の色分けが、着目SNPの３ディプロタイプに相当）。
１０万人、有病率１％、フェノコピーが２０％、遺伝率が６０％、着目座位は全遺伝成分の１％（ヒストグラムで、３ディプロタイプのヒストグラムにほとんど差がないことが、「全遺伝成分のわずかしかになっていないこと」と符合している）、リスクアレル頻度は１０％、HWE(F=0)の場合で、リスクアレルホモ（homo1)の有病者割合が0.016、ヘテロのそれが0.014、非リスクアレルホモ(homo2)のそれが、0.009となって、homo1/homo2のRRが1.79、hetero/homo2のRRが1.55になっていることがわかる。
単座位で考え、他の座位の影響を考えないとき、ホモ・ヘテロ・逆ホモに線形のリスクをモデルしているのに、homo1/homo2のRRとhetero/homo2のRRとがかなり近くなることもわかる
varg は着目SNPを除く遺伝因子の分散、vareは環境因子の分散、varmは着目SNPを除く遺伝因子と環境因子の和の分散、varsnpは着目SNPの分散、varaは着目SNPを含めた全遺伝因子と環境因子の和の分散

LocusRRSim<-function(N=10000,prevalence,phenocopy=0,heritability,locusFrac,alleleFreq,F=0,k=0){

#N<-10000000
#prevalence<-0.01
#phenocopy<-0.1
#print(c("N",N))
prevalence
#print(c("prevalence",prevalence))

phenocopy
#print(c("phenocopy",phenocopy))

heritability
#print(c("heritability",heritability))

locusFrac
#print(c("locusFrac",locusFrac))

alleleFreq
#print(c("alleleFreq",alleleFreq))

F
#print(c("F",F))


quantile<-1-prevalence*(1-phenocopy)
randphenocopy<-runif(N)
copycase<-which(randphenocopy<phenocopy*prevalence)

#gvfrac<-0.8
gvfrac<-heritability*(1-locusFrac)
#locusFrac<-0.05
#p<-0.2
p<-alleleFreq
g<-rnorm(N,mean=0,sd=sqrt(gvfrac))

#print(c("varg",var(g)))
e<-rnorm(N,mean=0,sd=sqrt(1-heritability))
#print(c("vare",var(e)))
m<-g+e
#print(c("varm",var(m)))
rand<-runif(N)
snprisk<-rep(0,N)

homo1<-which(rand<p*p+F*p*(1-p))
#hetero<-which((rand>=p*p)&&(rand<p*p+2*p*(1-p)))
homo2<-which(rand>=p*p+2*p*(1-p)-F*p*(1-p))

hetero<-setdiff(c(1:length(g)),homo1)
hetero<-setdiff(hetero,homo2)

X<-sqrt(heritability*locusFrac/(2*p*(1-p)))
snprisk[homo1]<-X*(1-p)*2
snprisk[homo2]<-X*(-p)*2
snprisk[hetero]<-X*(1-2*p+k)
#print(c("varsnp",var(snprisk)))
a<-m+snprisk
#print(c("vara",var(a)))

threshold<-quantile(a,quantile)
h<-hist(a)
hist(a[homo2],breaks=h$breaks,ylim=c(0,max(h$counts)),col="blue",density=30,main="",xlab="",ylab="")

par(new=T)
hist(a[hetero],breaks=h$breaks,ylim=c(0,max(h$counts)),col="purple",density=40,main="",xlab="",ylab="")
par(new=T)
hist(a[homo1],breaks=h$breaks,ylim=c(0,max(h$counts)),col="red",density=50,main="MAIN",xlab="XLAB",ylab="YLAB")

casehomo1<-length(union(which(a[homo1]>threshold),intersect(homo1,copycase)))
casehetero<-length(union(which(a[hetero]>threshold),intersect(hetero,copycase)))
casehomo2<-length(union(which(a[homo2]>threshold),intersect(homo2,copycase)))
conthomo1<-length(a[homo1])-casehomo1
conthetero<-length(a[hetero])-casehetero
conthomo2<-length(a[homo2])-casehomo2


casehomo1
#print(c("casehomo1",casehomo1))
#print(c("casehetero",casehetero))
#print(c("casehomo1",casehomo1))
#print(c("conthomo1",conthomo1))
#print(c("conthetero",conthetero))
#print(c("conthomo1",conthomo1))

casehetero
casehomo2
conthomo1
conthetero
conthomo2

prevhomo1<-casehomo1/(casehomo1+conthomo1)
prevhetero<-casehetero/(casehetero+conthetero)
prevhomo2<-casehomo2/(casehomo2+conthomo2)
prevhomo1
prevhetero
prevhomo2
#print(c("prevhomo1",prevhomo1))
#print(c("prevhetero",prevhetero))
#print(c("prevhomo1",prevhomo1))


RRhomo1<-prevhomo1/prevhomo2
RRhetero<-prevhetero/prevhomo2
#print(c("RRhomo1vshomo2",RRhomo1))
#print(c("RRheterovshomo2",RRhetero))

res<-list(N=N,prevalence=prevalence,phenocopy=phenocopy,heritability<-heritability,locusFrac=locusFrac,alleleFreq=alleleFreq,F=F,varg=var(g),vare=var(e),varm=var(m),varsnp=var(snprisk),vara=var(a),casehomo1=casehomo1,casehetero=casehetero,casehomo2=casehomo2,conthomo1=conthomo1,conthetero=conthetero,conthomo2=conthomo2,prevhomo1=prevhomo1,prevhetero=prevhetero,prevhomo2=prevhomo2,RRhomo1vshomo2=RRhomo1,RRheterovshomo2=RRhetero)
print(res)
return(res)


}



---------
N<-100000
prevalence<-0.01
phenocopy<-0.2
heritability<-0.6
locusFrac<-0.01
alleleFreq<-0.1
F<-0.0

res<-LocusRRSim(N=N,prevalence=prevalence,phenocopy=phenocopy,heritability=heritability,locusFrac=locusFrac,alleleFreq=alleleFreq,F=F)

$N
[1] 1e+05

$prevalence
[1] 0.01

$phenocopy
[1] 0.2

$heritability
[1] 0.6

$locusFrac
[1] 0.01

$alleleFreq
[1] 0.1

$F
[1] 0

$varg
[1] 0.5969842

$vare
[1] 0.4023782

$varm
[1] 0.9966114

$varsnp
[1] 0.005973443

$vara
[1] 1.001756

$casehomo1
[1] 22

$casehetero
[1] 228

$casehomo2
[1] 754

$conthomo1
[1] 1010

$conthetero
[1] 17770

$conthomo2
[1] 80216

$prevhomo1
[1] 0.02131783

$prevhetero
[1] 0.01266807

$prevhomo2
[1] 0.00931209

$RRhomo1vshomo2
[1] 2.289263

$RRheterovshomo2
[1] 1.36039



LocusRRSim(prevalence=prevalence,heritability=heritability,locusFrac=locusFrac,alleleFreq=alleleFreq,F=F)