- 分離比解析をやっている
- たくさんの家系データを作りたい
- 昨日の記事でグラフの成長ルールで家系を作る話を『なんちゃって家系図』として書いた
- 遺伝リスクやリスクアレルの伝達も含めて、少しきちんと作ってみる
- 設定
- 家系データの作り方
- 家系図を以下の要領で作成する
- trio行列は１行に１人
- 第１列は個人ID
- 第２列は父親ID(0は情報なし)
- 第３列は母親ID
- 第４列は性別(1:male,-1:female)
- 第５列は世代番号
- 手順1 第0世代にn.init人を作成する
- 手順2 初期n.init人の性別はランダムに与える
- 手順3 次いで、指定人数 n.iter人を次の要領で付け加える
- 手順4 すでに登録された個人から一人xを選ぶ
- 手順5-1 xが誰の親でもないときには、新たな配偶者を同世代に1人作成し、そのペアの間に、平均子供数mean.kidでポアソン分布に従う子供を作る
- 手順5-2すでに子供が居る場合で、親が居ない場合には、両親を１世代上に作成する
- 手順5-3 すでに子供が居て親もいるときは何もしない
- リスクローカスジェノタイプの与え方
- 一般集団のジェノタイプ頻度とジェノタイプ別リスクを与える
- 個人のリスクは全ローカスの単純和とする
- 環境要因を加えた上での有病者の指定の仕方
- 全個人のジェノタイプリスクの分散と指定した遺伝率とから環境要因の分散を逆算し、それを分散とする正規乱数を環境リスクとして加えて、個人の全リスクとする
- その上で総人数のうちの指定割合(有病率)を全リスクの上位から有病者として定める
- 有病者の選定には、もう少し工夫して、年齢を考慮したり、色々工夫してもよい
- 実際の利用にあたっては、有病者を見出し、そこから家系をたどる形で、発端者からの家系サンプリングなどをすることが可能。その際は、グラフオブジェクトに入れて、グラフ探索のアルゴリズムを使うのが便利なはず
my.pedigree <- function(n.init=3,n.iter=20,mean.kid=2){
trio <- matrix(0,n.init,5)
m <- length(trio[1,])
trio[1:n.init,1] <- 1:n.init
trio[1:n.init,4] <- sample(c(-1,1),n.init,replace=TRUE)
trio[,5] <- 0
cnt <- n.init
tmp.f <- function(m,k){
(m/(1-exp(-m))-k)^2
}
xmin <- optimize(tmp.f, c(0, mean.kid), tol = 0.0001, k = mean.kid)
mean.kid. <- xmin[[1]]
for(i in 1:n.iter){
s <- sample(trio[,1],1)
if(length(which(trio[,2:3]==s))==0){
trio <- rbind(trio,c(cnt+1,rep(0,m-1)))
dp <- dpois(1:100,mean.kid.)
dp <- dp/sum(dp)
num.kids <- sample(1:100,1,prob=dp)
tmp <- matrix(0,num.kids,m)
tmp[,1] <- cnt+1+(1:num.kids)
if(trio[s,4]==1){
tmp[,2] <- s
tmp[,3] <- cnt+1
}else{
tmp[,2] <- cnt+1
tmp[,3] <- s
}
tmp[,4] <- sample(c(-1,1),num.kids,replace=TRUE)
tmp[,5] <- trio[s,5] + 1
trio <- rbind(trio,tmp)
trio[cnt+1,4] <- -trio[s,4]
trio[cnt+1,5] <- trio[s,5]
cnt <- cnt+1+num.kids
}else if(trio[s,2] == 0 & trio[s,3] == 0){
trio <- rbind(trio,c(cnt+1,rep(0,m-1)))
trio <- rbind(trio,c(cnt+2,rep(0,m-1)))
trio[cnt+1,4] <- 1
trio[cnt+2,4] <- -1
trio[c(cnt+1,cnt+2),5] <- trio[s,5]-1
trio[s,2] <- cnt+1
trio[s,3] <- cnt+2
cnt <- cnt+2
}
}
trio
}
- 次に１ローカスのジェノタイプを家系図に合わせて発生する関数
make.genotype <- function(trio,prob){
ret <- rep(0,length(trio[,1]))
ord <- order(trio[,5])
for(i in 1:length(ord)){
if(trio[ord[i],2]==0 & trio[ord[i],3]==0){
ret[ord[i]] <- sample(c(0,1,2),1,prob=prob)
}else{
parents <- trio[ord[i],2:3]
pat <- ret[trio[ord[i],2]]/2
mat <- ret[trio[ord[i],3]]/2
if(pat==0.5){
pat <- sample(c(0,1),1)
}
if(mat==0.5){
mat <- sample(c(0,1),1)
}
ret[ord[i]] <- pat+mat
}
}
ret
}
- ついで、家系図情報、リスク座位のジェノタイプ頻度情報、リスク座位のジェノタイプ別リスク情報、遺伝率、有病率を与えて、個人のジェノタイプ行列、個人のジェノタイプリスク行列、個人の総リスクベクトル、疾病フェノタイプベクトルを返す関数
make.affected.pedigree <- function(trio,p.g,r.g,h,prev){
genotype <- matrix(0,length(trio[,1]),n.loci)
for(i in 1:n.loci){
genotype[,i] <- make.genotype(trio,p.g[i,])
}
R.mat <- genotype
for(i in 1:n.loci){
tmp <- r.g[i,]
R.mat[,i] <- tmp[genotype[,i]+1]
}
R.sum <- apply(R.mat,1,sum)
var.R <- var(R.sum)
var.E <- var.R*(1-herit)/herit
E <- rnorm(length(trio[,1]),0,sqrt(var.E))
R.E <- R.sum + E
affected <- rep(0,length(trio[,1]))
affected[which(R.E> quantile(R.E,1-prev))] <- 1
return(list(pedigree=trio,genotype=genotype,gen.risk.mat=R.mat,gen.risk=R.sum,affected=affected))
}
n.loci <- 10
p.a <- runif(n.loci)*0.1
p.g <- cbind(p.a^2,2*p.a*(1-p.a),(1-p.a)^2)
r.a <- c(10,rep(1,n.loci-1))
r.g <- cbind(r.a*2+1,r.a+1,rep(1,n.loci))
herit <- 0.5
prev <- 0.1
n.init <- 10
n.iter <- 50
mean.kid <-2
my.ped <- my.pedigree(n.init=n.init,n.iter=n.iter,mean.kid=mean.kid)
aff.pedigree <- make.affected.pedigree(my.ped,p.g,r.g,herit,prev)
library(kinship2)
plot.affected.pedigree <- function(ap){
sex <- ap$pedigree[,4]
sex[which(sex==-1)] <- 2
my.ped <- pedigree(ap$pedigree[,1],ap$pedigree[,2],ap$pedigree[,3],sex,ap$affected)
plot(my.ped,symbolosize=0.001)
}
- グラフオブジェクト化するために、エッジリストを作ればよい
el.from.trio <- function(trio){
non.zero <- which(apply((trio[,2:3])^2,1,sum)!=0)
el <- rbind(cbind(trio[non.zero,2],trio[non.zero,1]),cbind(trio[non.zero,3],trio[non.zero,1]))
el
}
library(igraph)
el <- el.from.trio(my.ped)
g <- graph.edgelist(el)
plot(g,vertex.color=aff.pedigree$pedigree[,4]+3)
plot(g,vertex.color=aff.pedigree$affected+1)
