少し真面目に家系図と家系データ

  • 分離比解析をやっている
  • たくさんの家系データを作りたい
  • 昨日の記事でグラフの成長ルールで家系を作る話を『なんちゃって家系図』として書いた
  • 遺伝リスクやリスクアレルの伝達も含めて、少しきちんと作ってみる
  • 設定
    • 家系データの作り方
      • 家系図を以下の要領で作成する
      • trio行列は1行に1人
      • 第1列は個人ID
      • 第2列は父親ID(0は情報なし)
      • 第3列は母親ID
      • 第4列は性別(1:male,-1:female)
      • 第5列は世代番号
      • 手順1 第0世代にn.init人を作成する
      • 手順2 初期n.init人の性別はランダムに与える
      • 手順3 次いで、指定人数 n.iter人を次の要領で付け加える
      • 手順4 すでに登録された個人から一人xを選ぶ
      • 手順5-1 xが誰の親でもないときには、新たな配偶者を同世代に1人作成し、そのペアの間に、平均子供数mean.kidでポアソン分布に従う子供を作る
      • 手順5-2すでに子供が居る場合で、親が居ない場合には、両親を1世代上に作成する
      • 手順5-3 すでに子供が居て親もいるときは何もしない
    • リスクローカスジェノタイプの与え方
      • 一般集団のジェノタイプ頻度とジェノタイプ別リスクを与える
      • 個人のリスクは全ローカスの単純和とする
    • 環境要因を加えた上での有病者の指定の仕方
      • 全個人のジェノタイプリスクの分散と指定した遺伝率とから環境要因の分散を逆算し、それを分散とする正規乱数を環境リスクとして加えて、個人の全リスクとする
      • その上で総人数のうちの指定割合(有病率)を全リスクの上位から有病者として定める
      • 有病者の選定には、もう少し工夫して、年齢を考慮したり、色々工夫してもよい
    • 実際の利用にあたっては、有病者を見出し、そこから家系をたどる形で、発端者からの家系サンプリングなどをすることが可能。その際は、グラフオブジェクトに入れて、グラフ探索のアルゴリズムを使うのが便利なはず

# trio行列は1行に1人
# 第1列は個人ID
# 第2列は父親ID(0は情報なし)
# 第3列は母親ID
# 第4列は性別(1:male,-1:female)
# 第5列は世代番号
# 第6列は情報(たとえば、着目アレルの本数)

# 1 第0世代にn.init人を作成する
# 初期n.init人の性別はランダムに与える
# 次いで、指定人数 n.iter人を次の要領で付け加える
# すでに登録された個人から一人xを選ぶ
# xが誰の親でもないときには、新たな配偶者を同世代に1人作成し
# そのペアの間に、平均子供数mean.kidでポアソン分布に従う子供を作る
# すでに子供が居る場合で、親が居ない場合には、両親を1世代上に作成する
# すでに子供が居て親もいるときは何もしない

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
	# 「子がいる場合」の「子の数の平均」がmean.kidなので、ポアソン分布に渡すパラメタはmean.kidより小さい
	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
}
  • 次に1ローカスのジェノタイプを家系図に合わせて発生する関数
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){
	# n.lociのジェノタイプを作る
	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))
}
  • さて、実行
    • 条件設定(1)
      • リスク座位のこと
# n.loci個のリスク座位
n.loci <- 10
# それらのアレル頻度とHWEを仮定したジェノタイプ頻度
## 適当にアレル頻度を発生させて
p.a <- runif(n.loci)*0.1
## HWEに従うジェノタイプ頻度を作成
p.g <- cbind(p.a^2,2*p.a*(1-p.a),(1-p.a)^2)

# 各座位のジェノタイプ別リスク
## 1強、n.loci-1弱のリスク設定
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
    • 条件設定(2)
      • 家系の作成
# 0世代の独立個人数
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)
  • プロットしてみる
# 性別が-1,1になっていることを修正しつつ
# kinship2パッケージの家系図描図関数を使う
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)