■

• こちらで複合遺伝性疾患の家系図を描くソースを作った
• 今回は、単一遺伝子病のそれを作る
• 常染色体劣性・優性
• Ｘ染色体劣性・優性
• Ｙ染色体
• ミトコンドリア

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

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

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
}
# my.pedigree()の出力trrioとアレル頻度を与えて、常染色体のジェノタイプを作る(ＨＷＥ仮定)
make.genotype.auto <- function(trio,p.a){
	prob <- c(p.a^2,2*p.a*(1-p.a),(1-p.a)^2)
	make.genotype(trio,prob)
}
# 同、Ｘ染色体ジェノタイプ
make.genotype.x <- function(trio,p.a){
	prob.male <- c(p.a,0,1-p.a)
	prob.female <- c(p.a^2,2*p.a*(1-p.a),(1-p.a)^2)
	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){
			if(trio[ord[i],4]==1){
				ret[ord[i]] <- sample(c(0,1,2),1,prob=prob.male)
			}else if(trio[ord[i],4]==-1){
				ret[ord[i]] <- sample(c(0,1,2),1,prob=prob.female)
			}
		}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)
			}
			if(trio[ord[i],4]==1){# male offspring
				ret[ord[i]] <- mat * 2
			}else if(trio[ord[i],4]==-1){# female offspring
				ret[ord[i]] <- mat + pat
			}
		}
	}
	ret
}
# 同、Ｙ染色体ジェノタイプ
make.genotype.y <- function(trio,p.a){
	prob.male <- c(p.a,0,1-p.a)
	prob.female <- c(1,0,0)
	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){
			if(trio[ord[i],4]==1){
				ret[ord[i]] <- sample(c(0,1,2),1,prob=prob.male)
			}else if(trio[ord[i],4]==-1){
				ret[ord[i]] <- sample(c(0,1,2),1,prob=prob.female)
			}

		}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)
			#}
			if(trio[ord[i],4]==1){# male offspring
				ret[ord[i]] <- pat * 2
			}else if(trio[ord[i],4]==-1){# female offspring
				#ret[ord[i]] <- mat + pat
				ret[ord[i]] <- 0
			}
		}
	}
	ret
}
# 同、ミトコンドリアジェノタイプ(ただし、ミトコンドリアは１人に１タイプ
make.genotype.mit <- function(trio,p.a){
	prob.male <- c(p.a,0,1-p.a)
	prob.female <- c(p.a,0,1-p.a)
	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){
			if(trio[ord[i],4]==1){
				ret[ord[i]] <- sample(c(0,1,2),1,prob=prob.male)
			}else if(trio[ord[i],4]==-1){
				ret[ord[i]] <- sample(c(0,1,2),1,prob=prob.female)
			}

		}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)
			#}
			if(trio[ord[i],4]==1){# male offspring
				ret[ord[i]] <- mat * 2
			}else if(trio[ord[i],4]==-1){# female offspring
				ret[ord[i]] <- mat * 2
			}
		}
	}
	ret
}
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
}
# ジェノタイプ値は0,1,2として、それに重みをつけてスコアリングする
my.phenotype.score <- function(g,w=c(0,0.5,1)){
	ret <- rep(0,length(g))
	for(i in 1:3){
		ret[g==(i-1)] <- ret[g==(i-1)] + w[i]
	}
	ret
}
# 優性モデルでのスコアリング(常染色体もＸ染色体も使える。Ｙとミトコンドリアにも使える)
# comple=TRUEは、表現型の「裏」を描図
my.phenotype.score.dom <- function(g,comple=FALSE){
	if(comple){
		return(my.phenotype.score(g,w=c(1,0,0)))
	}
	my.phenotype.score(g,w=c(0,1,1))
}
# 劣性モデルでのスコアリング
my.phenotype.score.rec <- function(g,comple=FALSE){
	if(comple){
		return(my.phenotype.score(g,w=c(1,1,0)))
	}
	my.phenotype.score(g,w=c(0,0,1))
}

# my.pedigree()は男女を1,-1とするが、kinship2パッケージは1,2とするので、その変換
my.change.sex.val <- function(s){
	s[which(s==-1)] <- 2
	return(s)
}
# kinshipt2パッケージのpedigree()関数とそのplot.pedigree()関数を使う
library(kinship2)
my.plot.pedigree <- function(my.ped,af,sex.change=TRUE,s.size=1){
	if(sex.change){
		my.ped[,4] <- my.change.sex.val(my.ped[,4])
	}
	my.ped <- pedigree(my.ped[,1],my.ped[,2],my.ped[,3],my.ped[,4],af)
	plot(my.ped,symbolsize=s.size)	
}

# 0世代の独立個人数
n.init <- 5
# 家系図を膨らませる処理ステップ数
n.iter <- 20
# 平均子供数
mean.kid <- 5

my.ped <- my.pedigree(n.init=n.init,n.iter=n.iter,mean.kid=mean.kid)

af.auto <- 0.8

g.auto <- make.genotype.auto(my.ped,af.auto)

par(mfcol=c(1,2))
my.plot.pedigree(my.ped,af=my.phenotype.score.dom(g.auto))
my.plot.pedigree(my.ped,af=my.phenotype.score.dom(g.auto,comple=TRUE))
af.x <- 0.8

g.x <- make.genotype.x(my.ped,af.x)

my.plot.pedigree(my.ped,af=my.phenotype.score.dom(g.x))
my.plot.pedigree(my.ped,af=my.phenotype.score.dom(g.x,comple=TRUE))
my.plot.pedigree(my.ped,af=my.phenotype.score.rec(g.x))
my.plot.pedigree(my.ped,af=my.phenotype.score.rec(g.x,comple=TRUE))


af.y <- 0.3


g.y <- make.genotype.y(my.ped,af.y)
my.plot.pedigree(my.ped,af=my.phenotype.score.dom(g.y))
my.plot.pedigree(my.ped,af=my.phenotype.score.dom(g.y,comple=TRUE))

af.mit <- 0.8
g.mit <- make.genotype.mit(my.ped,af.mit)
my.plot.pedigree(my.ped,af=my.phenotype.score.dom(g.mit))
my.plot.pedigree(my.ped,af=my.phenotype.score.dom(g.mit,comple=TRUE))