DNA二重らせんの高次畳みこみ - ryamadaの遺伝学・遺伝統計学メモ

昨日の記事で高次らせんを描いた
これを一昨日の記事のDNA２重らせんと組み合わせると…
library(rgl)
# 基本的には高次らせん
ori <- c(0,0,0)
# 第i円の半径。DNA２重らせんの円周は10オングストロームに合わせ、Rsの最終要素の値は10
Rs <- c(640,160,40,10)
# DNA２重らせんの１周あたりの伸びが34オングストロームなので、後ろから２番目の円の半径x角に関してはそれを反映させた値とするべくTs[length(Rs)-1]をRs[length(Rs)-1]の値で調整している
# 後で出るパラメタtが１刻みの整数であるので、Tsの最終円の値は1\
Ts <- rep(1,length(Rs))
Ts[length(Rs)-3] <- 0.005
Ts[length(Rs)-2] <- 0.05
Ts[length(Rs)-1] <- 34/Rs[length(Rs)-1]
# 高次階層数
n.steps <- length(Rs)
max.t <- 0.1
# 最終円の１周あたりの塩基数を指定
n.base <- 10.4
t.n <- 1000
t <- (0:(t.n-1))/n.base*2*pi

unit.vector <- function(X){
	v <- sqrt(apply(X^2,1,sum))
	X/v
}
outer.prod <- function(v1,v2){
	c(v1[2]*v2[3]-v1[3]*v2[2],v1[3]*v2[1]-v1[1]*v2[3],v1[1]*v2[2]-v1[2]*v2[3])
}
outer.prod.V <- function(V1,V2){
	cbind(V1[,2]*V2[,3]-V1[,3]*V2[,2],V1[,3]*V2[,1]-V1[,1]*V2[,3],V1[,1]*V2[,2]-V1[,2]*V2[,1])
}
X <- list()
V <- list()

X[[1]] <- V[[1]] <- cbind(rep(0,t.n),rep(0,t.n),rep(0,t.n))
X[[2]] <- cbind(Rs[1]*cos(Ts[1]*t),Rs[1]*sin(Ts[1]*t),rep(0,t.n))
V[[2]] <- cbind(-Rs[1]*sin(Ts[1]*t),Rs[1]*cos(Ts[1]*t),rep(0,t.n))
comple.X <- matrix(0,t.n,3)
for(i in 2:n.steps){
	tmp.x <- X[[i]]-X[[i-1]]
	tmp.x <- unit.vector(tmp.x)
	tmp.v <- unit.vector(V[[i]])
	tmp.z <- outer.prod.V(tmp.x,tmp.v)
	tmp.X <- cbind(Rs[i]*cos(Ts[i]*t),Rs[i]*sin(Ts[i]*t),rep(0,t.n))
	V[[i+1]] <- cbind(-Rs[i]*sin(Ts[i]*t),Rs[i]*cos(Ts[i]*t),rep(0,t.n))
	for(j in 1:t.n){
		tmp.X[j,] <- cbind(tmp.x[j,],tmp.z[j,],tmp.v[j,]) %*% tmp.X[j,]
		V[[i+1]][j,] <- cbind(tmp.x[j,],tmp.z[j,],tmp.v[j,]) %*% V[[i+1]][j,]
	}
	X[[i+1]] <- X[[i]] + tmp.X
	if(i==n.steps){
		tmp.X <- cbind(Rs[i]*cos(Ts[i]*t+pi),Rs[i]*sin(Ts[i]*t+pi),rep(0,t.n))
		for(j in 1:t.n){
			tmp.X[j,] <- cbind(tmp.x[j,],tmp.z[j,],tmp.v[j,]) %*% tmp.X[j,]
		}
		comple.X <- X[[i]] + tmp.X
	}
}

plot3d.minmax <- function(X){
	min.X <- min(X)
	max.X <- max(X)
	tmp <- rbind(X,rep(min.X,length(X[1,])),rep(max.X,length(X[1,])))
	plot3d(tmp,col="white",axes=FALSE)
}

means <- apply(comple.X,2,mean)
comple.X <- t(t(comple.X)-means)
X[[n.steps+1]] <- t(t(X[[n.steps+1]])-means)

tmp.d <- dist(X[[n.steps+1]])
max.X <- max(X[[n.steps+1]])
min.X <- min(X[[n.steps+1]])
dd <- max.X-min.X


col <- sample(0:3,length(t),replace=TRUE)
new.col <- (col+2)%%4
# 塩基に彩色した球を対応づける
base.col <- c("red","green","blue","yellow")

plot3d.minmax(X[[n.steps+1]])
#spheres3d(X[[n.steps+1]],col=base.col[col+1],radius=max(min(tmp.d)*0.5,dd*0.001))
spheres3d(X[[n.steps+1]],col=base.col[col+1],radius=2)
s <- rep(1:t.n,each=2)
s <- s[2:(length(s)-1)]
segments3d(X[[n.steps+1]][s,],col=gray(0.1))

#spheres3d(comple.X,col=base.col[new.col+1],radius=max(min(tmp.d)*0.5,dd*0.001))
spheres3d(comple.X,col=base.col[new.col+1],radius=2)
segments3d(comple.X[s,],col=gray(0.1))

tandem.xyz <- rbind(X[[n.steps+1]],comple.X)
bridges <- c(t(matrix(1:(2*t.n),ncol=2)))
segments3d(tandem.xyz[bridges,],col=gray(0.7))