染色体グラフ

  • こちらで、ランダムグラフを勉強している
  • Y染色体の伝搬をランダムグラフで扱ったり、巨大な連結グラフに確率過程を持ち込むことで、数学的に扱って、グラフの特徴を定性的・定量的に求める考え方が書かれている
  • 染色体の伝搬をY染色体のように単純に考えずに組換えありで扱うと、どんな風になるか(扱えるのかも含めて)、夏休みらしく、うまく行かなくても、いいか、というスタンスでやってみることにする
  • 変数設定
    • グラフの構成成分
      • ノードは、すべての塩基とする
      • エッジは、染色体上のつながりと、親子関係の伝達との2種類とする
    • ノードの属性
      • ノードは座標を持つ
        • 染色体をこの世という空間にある点とみなしたときの、その位置座標
        • 染色体上の塩基番号という座標
          • 染色体には、必ず交叉が起きるというルールがあるので、染色体に必ず交叉を入れるような設定にして、染色体間の距離は、そのような交叉発生のルールのもとで、ほぼ0.5の確率で奇数回の交叉が入るようなものを取ればよさそうだ
        • 時刻という座標
    • エッジの属性
      • 染色体というつながりか、伝達というつながりか、の2種類に分かれる
  • そうは言っても
    • まず、染色体のグラフを作って、その上で、塩基のノードとそれらの関連にエッジをつけるのが手順としてはよさそうだ
  • まずはイザナギイザナミ(アダムとイブ)からスタートにしよう
    • 世代ごとに少しずつ人口をふやしつつ
    • ペアリングに制約を入れることを考慮し
    • 移住しながら…

# 第1世代は2人
init.n <- 2
# 世代数
n.gen <- 50
# 世代別人口
pop.gen <- c(init.n)
# 空間座標
loc <- list()
loc[[1]] <- matrix(0,pop.gen[1],2)
loc[[1]] <- loc[[1]] + rnorm(length(loc[[1]]))
# 親情報
parents <- list()
parents[[1]] <- matrix(0,pop.gen[1],2)
# 空間的に近すぎる個人ペアはペアリング不可
thres = 0 # 0とは「同一個人を両親とするのは不可」の意味
# 世代ごとの人口増加率
incr.ratio <- 1.05
# 移住パラメタ
migration <- runif(2)

for(i in 2:n.gen){
	# ペアリングは近い方がペアになり易い
	d <- matrix(0,pop.gen[i-1],pop.gen[i-1])
	for(j in 1:(pop.gen[i-1]-1)){
		for(j2 in (j+1):pop.gen[i-1]){
			d[j,j2] <- d[j2,j] <- sum((loc[[i-1]][j,]-loc[[i-1]][j2,])^2)
		}
	}
	d[which(d <= thres)] <- Inf
	prob <- exp(-d^2/2)
	prob <- prob/sum(prob)
	# 次世代人口の最大値
	max.child <- pop.gen[i-1]*incr.ratio * 5
	# 子供数は、2から最大値まで(絶滅しない設定)
	n.child <- sample(2:max.child,1,prob = dpois(2:max.child,pop.gen[i-1]*incr.ratio))
	pop.gen <- c(pop.gen, n.child)
	#print(pop.gen)
	parents.number <- sample(1:length(d),pop.gen[i],replace=TRUE,prob = prob)
	address <- matrix(1:length(d),pop.gen[i-1],pop.gen[i-1])
	parents.list <- which(address <= max(address),arr.ind=TRUE)
	parents[[i]] <- parents.list[parents.number,]
	#print(parents[[i]])
	loc[[i]] <- (loc[[i-1]][parents[[i]][,1],] + loc[[i-1]][parents[[i]][,2],])/2
	loc[[i]] <- loc[[i]] + cbind(rnorm(length(loc[[i]][,1]),rnorm(1,migration[1])),rnorm(length(loc[[i]][,2]),rnorm(1,migration[2])))
}

xlim <- ylim <- range(unlist(loc))
plot(loc[[1]],xlim = xlim, ylim = ylim,col=1,cex = 0.5,pch =20)
for(i in 2:n.gen){
	#plot(loc[[i]],xlim = xlim,ylim=ylim)
	points(loc[[i]],col=i,cex = 0.5,pch=20)
	Sys.sleep(0.1)
}
  • 少し遊んでみる
library(MCMCpack)
# 第1世代は2人
init.n <- 10
# 世代数
n.gen <- 100
# 世代別人口
pop.gen <- c(init.n)
# 空間座標
loc <- list()
loc[[1]] <- matrix(0,pop.gen[1],2)
loc[[1]] <- loc[[1]] + rnorm(length(loc[[1]]),sd=1)
#loc[[1]] <- loc[[1]] + rdirichlet(init.n,rep(1,2))*10
#loc[[1]] <- loc[[1]] - 0.5
# 親情報
parents <- list()
parents[[1]] <- matrix(0,pop.gen[1],2)
# 空間的に近すぎる個人ペアはペアリング不可
#thres = 0 # 0とは「同一個人を両親とするのは不可」の意味
thres = mean(abs((c(loc[[1]])))) # 0とは「同一個人を両親とするのは不可」の意味
# 世代ごとの人口増加率
incr.ratio <- 1.01
# 移住パラメタ
migration <- runif(2)

for(i in 2:n.gen){
	# ペアリングは近い方がペアになり易い
	d <- matrix(0,pop.gen[i-1],pop.gen[i-1])
	for(j in 1:(pop.gen[i-1]-1)){
		for(j2 in (j+1):pop.gen[i-1]){
			d[j,j2] <- d[j2,j] <- sum((loc[[i-1]][j,]-loc[[i-1]][j2,])^2)
		}
	}
	print(length(which(d <= thres)))
	d[which(d <= thres)] <- 1/d[which(d <= thres)]
	prob <- exp(-d^2/2)
	prob <- prob/sum(prob)
	# 次世代人口の最大値
	max.child <- pop.gen[i-1]*incr.ratio * 5
	# 子供数は、2から最大値まで(絶滅しない設定)
	n.child <- sample(2:max.child,1,prob = dpois(2:max.child,pop.gen[i-1]*incr.ratio))
	pop.gen <- c(pop.gen, n.child)
	#print(pop.gen)
	parents.number <- sample(1:length(d),pop.gen[i],replace=TRUE,prob = prob)
	address <- matrix(1:length(d),pop.gen[i-1],pop.gen[i-1])
	parents.list <- which(address <= max(address),arr.ind=TRUE)
	parents[[i]] <- parents.list[parents.number,]
	#print(parents[[i]])
	loc[[i]] <- (loc[[i-1]][parents[[i]][,1],] + loc[[i-1]][parents[[i]][,2],])/2
	loc[[i]] <- loc[[i]] + cbind(sign(loc[[i]][,1])*rnorm(length(loc[[i]][,1]),rnorm(1,migration[1])),sign(loc[[i]][,2])*rnorm(length(loc[[i]][,2]),rnorm(1,migration[2])))
	loc[[i]] <- loc[[i]] + cbind((loc[[i]][,2]),-(loc[[i]][,1])) /mean(abs(loc[[i]]))
}

xlim <- ylim <- range(unlist(loc))
plot(loc[[1]],xlim = xlim, ylim = ylim,col=1,cex = 1,pch =20)
for(i in 2:n.gen){
	#plot(loc[[i]],xlim = xlim,ylim=ylim)
	points(loc[[i]],col=i,cex = 1,pch=20)
	for(j in 1:pop.gen[i]){
		segments(loc[[i]][j,1],loc[[i]][j,2],loc[[i-1]][parents[[i]][j,1],1],loc[[i-1]][parents[[i]][j,1],2])
		segments(loc[[i]][j,1],loc[[i]][j,2],loc[[i-1]][parents[[i]][j,2],1],loc[[i-1]][parents[[i]][j,2],2])
		Sys.sleep(0.01)
	}
	Sys.sleep(0.1)
}