2013-12-24

DNA２重らせんを使ってknitrに慣れる

R knitr rgl

knitrを使ってR Markdown文書の扱いを練習している。
DNA２重らせんを描いていみよう

電子書籍化してみました

作者: ryamada
発売日: 2013/12/30
メディア: Kindle版
この商品を含むブログ (1件) を見る

library(rgl)
# 二重らせん
# 塩基対の数
n.pt <- 50
# らせん１周当たりの塩基数
n <- 10.4
# 円半径
r1 <- 10
# 円１周あたりのらせんの伸び
r2 <- 34
# 塩基ごとの円周の角
t <- (0:(n.pt-1))/n*2*pi
# 座標
x <- r1*cos(t)
y <- r1*sin(t)
z <- t/(2*pi)*r2
xyz <- cbind(x,y,z)
# 隣接塩基を線分で結ぶ
# 隣接する塩基同士を線分でつなぐための工夫
# (1,2,2,3,3,4,4,5,...)というベクトルを作る
s <- rep(1:n.pt,each=2)
s <- s[2:(length(s)-1)]
segments3d(xyz[s,])
# 塩基を0,1,2,3の４タイプで指定する。ランダム
col <- sample(0:3,length(t),replace=TRUE)
# 塩基に彩色した球を対応づける
base.col <- c("red","green","blue","yellow")
base.size <- 2
# 3次元プロット
xyz. <- rbind(xyz,rep(min(xyz),3),rep(max(xyz),3))
plot3d(xyz.,col=rep("white",length(xyz.[,1])),axes=FALSE)
spheres3d(xyz,col=base.col[col+1],radius=base.size)
segments3d(xyz[s,])
# 角piだけ回したらせんの点の座標とその隣接塩基間を結ぶ線分を描きこむ
new.x <- r1*cos(t+pi)
new.y <- r1*sin(t+pi)
new.z <- z
new.xyz <- cbind(new.x,new.y,new.z)
new.s <- s
new.col <- (col+2)%%4
spheres3d(new.xyz,col=base.col[new.col+1],radius=base.size)
segments3d(new.xyz[s,])

tandem.xyz <- rbind(xyz,new.xyz)
bridges <- c(t(matrix(1:(2*n.pt),ncol=2)))
segments3d(tandem.xyz[bridges,],col=gray(0.6))

# DNA二重らせんを使ってRの練習をする

DNA二重らせんの３次元プロットをしてみる。

## らせん

らせんとはどういうものだろうか。

２次元平面に円を描きながら、その２次元平面に垂直な方向に
等速で移動したときにできる軌跡がらせんである。

円をxy平面に描き、z軸方向に等速で動くこととする。

円の１周は$2\pi$なので、１周当たり、n個の点を打つとすれば、
\(
\theta_i = \frac{i \times 2\pi}{n}; i=0,1,2,...
\)

のような角を与えればよい。

円周の半径を$r_1$とすれば

\(
x_i = r_1 \times \cos(\theta_i)
\)

\(
y_i = r_1 \times \sin(\theta_i)
\)

DNA二重らせんの場合、円の半径は10オングストローム、
１周あたりの塩基数は約10.4と言われているので、その値を使う。

Rを使って描いてみる
```{r,fig.width=7,fig.height=7}
# すべての点の数をn.ptとする
n.pt <- 15
i.s <- 0:(n.pt-1)
# １周当たりの点の数
n <- 10.4
# 角
theta <- i.s * 2 * pi / n
# 円の半径
r1 <- 10
# x,y座標
x <- r1 * cos(theta)
y <- r1 * sin(theta)
# プロットしてみる
plot(x,y,type="b")
```

次にz軸に等速移動させる。

円を１周するごとに$r_2$進むとすると、
角$\theta$のときのz軸方向の移動距離は

\(
\frac{\theta}{2\pi}\times r_2
\)

で与えられる。

DNA二重らせんの場合$r_2=34$オングストロームと言われている。
Rを使って描いてみる。

```{r,fig.width=7,fig.height=7}
# らせんのz方向測度
r2 <- 34
z <- theta/(2*pi)*r2
plot(x,z,type="b")
abline(h=r2,col=2)
```

$z=0$から$z=r_2$まで進んだところでx座標に関してひとめぐりしている。

これを3次元プロットしてみる。
```{r setup, rgl = TRUE}
library(rgl)
knit_hooks$set(rgl = hook_rgl)
plot3d(x,y,z)
```

## 見やすい工夫をする

らせん上の点は描けたけれど、点と点の間を線分で結ぶことにする。

```{r, rgl=TRUE}
# 塩基対の数
n.pt <- 50
# らせん１周当たりの塩基数
n <- 10.4
# 円半径
r1 <- 10
# 円１周あたりのらせんの伸び
r2 <- 34
# 塩基ごとの円周の角
t <- (0:(n.pt-1))/n*2*pi
# 座標
x <- r1*cos(t)
y <- r1*sin(t)
z <- t/(2*pi)*r2
xyz <- cbind(x,y,z)

# 塩基を0,1,2,3の４タイプで指定する。ランダム
col <- sample(0:3,length(t),replace=TRUE)

# 3次元プロット
library(rgl)
knit_hooks$set(rgl = hook_rgl)
plot3d(xyz)
# 塩基に彩色した球を対応づける
base.col <- c("red","green","blue","yellow")
base.size <- 2
spheres3d(xyz,col=base.col[col+1],radius=base.size)
# 隣接塩基を線分で結ぶ
# 隣接する塩基同士を線分でつなぐための工夫
# (1,2,2,3,3,4,4,5,...)というベクトルを作る
s <- rep(1:n.pt,each=2)
s <- s[2:(length(s)-1)]
segments3d(xyz[s,])
```

さらにもう少し工夫をする。

x,y,z軸のスケールが異なっているので、これを同じにする。

そのために、以下では、xyzの最小値をx,y,z座標に持つ点と、xyzの最大値をx,y,z座標に持つ点とを描図オブジェクトに加え、それは「白色」で見えないように打点している。

また、座標の枠も消しておく。

```{r,rgl=TRUE}
knit_hooks$set(rgl = hook_rgl)
xyz. <- rbind(xyz,rep(min(xyz),3),rep(max(xyz),3))
plot3d(xyz.,col=rep("white",length(xyz.[,1])),axes=FALSE)
spheres3d(xyz,col=base.col[col+1],radius=base.size)
segments3d(xyz[s,])
```

## ２重らせん

DNAは２重らせんである。

DNA２重らせんでは２本のらせんが並行にパイプのような構造を作っている。

実際、座標を与えるときには、円周に関して角$\pi$だけ回転したらせんを付け加えるとできる。

付け加える鎖は相補鎖のことである。

Rでやってみる

```{r,rgl=TRUE}
knit_hooks$set(rgl = hook_rgl)
xyz. <- rbind(xyz,rep(min(xyz),3),rep(max(xyz),3))
plot3d(xyz.,col=rep("white",length(xyz.[,1])),axes=FALSE)
spheres3d(xyz,col=base.col[col+1],radius=base.size)
segments3d(xyz[s,])
# 角piだけ回したらせんの点の座標とその隣接塩基間を結ぶ線分を描きこむ
new.x <- r1*cos(t+pi)
new.y <- r1*sin(t+pi)
new.z <- z
new.xyz <- cbind(new.x,new.y,new.z)
new.s <- s
spheres3d(new.xyz,radius=base.size)
segments3d(new.xyz[s,])
```

## 相補結合と剰余

前の図では相補鎖の塩基の色を黒とした。

実際には、初めの鎖の塩基A,T,G,CにはT,A,C,Gが対応づく。

今、塩基を適当に決めるにあたって、0,1,2,3の数字をランダムに対応づけている。

これを用いる。

(A=0,G=1,T=2,C=3)とすれば、相補鎖に関しては、(0(A)->2(T),1(G)->3(C),2(T)->0(A),3(C)->1(G))という関係になる。

ここで、4で割った剰余のことを考えれば、
0+2を4で割った剰余は2、1+2を4で割った剰余は3、2+2を4で割った剰余は0、3+2を4で割った剰余は1であるから、
(0,1,2,3)に2を加えて、その4で割った剰余を取ることにすればｇ、相補塩基が決まる。

以下では、それを利用して、相補鎖の塩基の色を与えている。

```{r,rgl=TRUE}
knit_hooks$set(rgl = hook_rgl)
xyz. <- rbind(xyz,rep(min(xyz),3),rep(max(xyz),3))
plot3d(xyz.,col=rep("white",length(xyz.[,1])),axes=FALSE)
spheres3d(xyz,col=base.col[col+1],radius=base.size)
segments3d(xyz[s,])
# 角piだけ回したらせんの点の座標とその隣接塩基間を結ぶ線分を描きこむ
new.x <- r1*cos(t+pi)
new.y <- r1*sin(t+pi)
new.z <- z
new.xyz <- cbind(new.x,new.y,new.z)
new.s <- s
new.col <- (col+2)%%4
spheres3d(new.xyz,col=base.col[new.col+1],radius=base.size)
segments3d(new.xyz[s,])
```

これに、相補関係にある塩基間を結ぶ線分を加える。

元の鎖の塩基の座標行列と相補鎖の塩基の座標行列を連結し、
その行列の行番号の列(1,n.pt+1,2,n+pt+2,...)を用いて、線分を描きいれたのが以下。

```{r,rgl=TRUE}
knit_hooks$set(rgl = hook_rgl)
xyz. <- rbind(xyz,rep(min(xyz),3),rep(max(xyz),3))
plot3d(xyz.,col=rep("white",length(xyz.[,1])),axes=FALSE)
spheres3d(xyz,col=base.col[col+1],radius=base.size)
segments3d(xyz[s,])
# 角piだけ回したらせんの点の座標とその隣接塩基間を結ぶ線分を描きこむ
new.x <- r1*cos(t+pi)
new.y <- r1*sin(t+pi)
new.z <- z
new.xyz <- cbind(new.x,new.y,new.z)
new.s <- s
new.col <- (col+2)%%4
spheres3d(new.xyz,col=base.col[new.col+1],radius=base.size)
segments3d(new.xyz[s,])

tandem.xyz <- rbind(xyz,new.xyz)
bridges <- c(t(matrix(1:(2*n.pt),ncol=2)))
segments3d(tandem.xyz[bridges,],col=gray(0.6))
```
## まとめ

最後にコマンドをひとそろい書いておこう。

```{r}
# 二重らせん
# 塩基対の数
n.pt <- 50
# らせん１周当たりの塩基数
n <- 10.4
# 円半径
r1 <- 10
# 円１周あたりのらせんの伸び
r2 <- 34
# 塩基ごとの円周の角
t <- (0:(n.pt-1))/n*2*pi
# 座標
x <- r1*cos(t)
y <- r1*sin(t)
z <- t/(2*pi)*r2
xyz <- cbind(x,y,z)
# 隣接塩基を線分で結ぶ
# 隣接する塩基同士を線分でつなぐための工夫
# (1,2,2,3,3,4,4,5,...)というベクトルを作る
s <- rep(1:n.pt,each=2)
s <- s[2:(length(s)-1)]
segments3d(xyz[s,])
# 塩基を0,1,2,3の４タイプで指定する。ランダム
col <- sample(0:3,length(t),replace=TRUE)
# 塩基に彩色した球を対応づける
base.col <- c("red","green","blue","yellow")
base.size <- 2
# 3次元プロット
xyz. <- rbind(xyz,rep(min(xyz),3),rep(max(xyz),3))
library(rgl)
plot3d(xyz.,col=rep("white",length(xyz.[,1])),axes=FALSE)
spheres3d(xyz,col=base.col[col+1],radius=base.size)
segments3d(xyz[s,])
# 角piだけ回したらせんの点の座標とその隣接塩基間を結ぶ線分を描きこむ
new.x <- r1*cos(t+pi)
new.y <- r1*sin(t+pi)
new.z <- z
new.xyz <- cbind(new.x,new.y,new.z)
new.s <- s
new.col <- (col+2)%%4
spheres3d(new.xyz,col=base.col[new.col+1],radius=base.size)
segments3d(new.xyz[s,])

tandem.xyz <- rbind(xyz,new.xyz)
bridges <- c(t(matrix(1:(2*n.pt),ncol=2)))
segments3d(tandem.xyz[bridges,],col=gray(0.6))
```