塩基配列上での変異とアミノ酸置換

R 遺伝医学 数理モデル Rpackage Rを使いこなす

MikuHatsune2013-04-04

J Hepatol. 2009 May;50(5):895-905
HBVの薬剤応答具合と逆転写酵素領域の配列の変化をモデル化して考えた論文。
治療開始後48週間後では、反応が良かった群で治療開始前の塩基配列と比べて複雑度と多様度が減少している、ということらしい。
その数理モデルをシミュレーションしてみる。
　
Complexity
シャノンエントロピー
クローン数が全部でN個の配列で、配列s_iがp_iの割合で存在しているとすると
S_n=-\sum_i^n \frac{p_i log{p_i}}{log N}
と計算される。
　
Diversity
Hamming距離
ベクトルのうち、異なる要素があればそこの距離を1にしてしまうやり方。
配列にしてしまえば、異なる文字であれば1にして、和をとる。
vwrパッケージでできる。

library(vwr)
slen <- 20
set.seed(831)
s0 <- paste(sample(c("A", "G", "C", "T"), size=slen, replace=TRUE), sep="", collapse="")
s1 <- paste(sample(c("A", "G", "C", "T"), size=slen, replace=TRUE), sep="", collapse="")
s0; s1; hamming.distance(s0, s1)

[1] "CGTTACGATGTGATTAACCC"
[1] "GTCGAGGTAATAATTGGGTA"
GTCGAGGTAATAATTGGGTA
                  14

　
\frac{Ka}{Ks}比というモデルがある。
これは、塩基配列中のコドンでの塩基変異が、アミノ酸置換を起こす(non-synonymous, Ka or dN)か起こさないか(synonymous, Ks or dS)という塩基配列上でのポジションごとにいろいろ計算する、というもの。
ググっても解説してくれているサイトはあまりない気がする。検索するとけっこう上位に出てくる。
Rではseqinrパッケージで計算できる。

　
いま、適当な長さの配列を発生させたとする。
この配列の長さは、タンパクに翻訳されるように3の倍数にしてある。
この配列に、配列の長さのうち%の塩基がこれまた適当に塩基変異するとする。ただし、G→Gも塩基変化としてしまうことにする。
塩基変化する%をいろいろ変えてみて、KaとKsがどう変化するかをシミュレーションした。

library(seqinr)
# seqinr のexample
s <- read.alignment(file = system.file("sequences/test.phylip", package = "seqinr"), format = "phylip")
s0 <- kaks(s)
aa <- lapply(lapply(lapply(s$seq, s2c), translate), paste, sep="", collapse="") # s2c で文字列をバラバラにして、translate で翻訳する。
aa_aln <- as.alignment(5, s$nam, aa) # alignment でalignment オブジェクトにする。

AAs <- 100 # アミノ酸数
#HS <- s2c(s$seq[[3]]) 
HS <- sample(c("a", "c", "g", "t"), size=3*AAs, replace=TRUE) # 配列を適当に作成。
n1 <- 100 # 変異を起こす配列数
mutation_base <- length(HS) # 変異を起こす塩基数
niter <- 1 # シミュレーションしようと思っていたけど作成する配列数を増やせばいいことに気づいた。
res <- array(0, c(mutation_base, n1 - 1, 2)) #ka と Ks
ent <- rep(0, mutation_base) # シャノンエントロピー
for(m0 in seq(mutation_base)){
	mHS <- vector("list", n1 - 1)
	for(i in seq(mHS)){
		#m0 <- 2 #sample(length(HS), 1)
		mHS[[i]] <- replace(HS, sample(seq(HS), m0), sample(c("a", "g", "c", "t"), m0, replace=TRUE))
	}
	mHS <- lapply(c(list(HS), mHS), paste, sep="", collapse="")
	t0 <- table(unlist(mHS))
	ent[m0] <- -sum((t0/sum(t0))*log(t0/sum(t0))/sum(t0))
	s1 <- kaks(as.alignment(length(mHS), seq=mHS))
	res[m0, , 1] <- tail(as.matrix(s1$ka)[1, ], -1L)
	res[m0, , 2] <- tail(as.matrix(s1$ks)[1, ], -1L)
}
par(mfcol=c(1, 2))
for(i in seq(2)){
	boxplot(t(res[ , , i]), xaxt="n", xlab="proportion of mutated bases", ylab="value", ylim=c(0, 4), cex=0.5)
	title(c("Ka", "Ks")[i])
	axis(1, at=seq(0, 3*AAs, length=11), labels=seq(0, 1, length=11))
}

KaとKsはそれぞれ規格化されているので、適当な塩基変異の場合、\frac{Ka}{Ks}は1になる。
\frac{Ka}{Ks}>1の場合、タンパク質の変化を伴うような変異、つまりnon-synonymousな変異が積極的に起こっていると考えられ、生き残りをかけているような正の選択がなされていると考える。これは、未知の細菌やウイルスといった抗原に対する免疫反応を多様化させるといったところで役に立つ。
\frac{Ka}{Ks}<1の場合はその逆で、タンパク質変化を伴わないsynonymousな変異で保守的なんだろう。

図は配列中の変異する塩基数が多くなると、Ka, Ksともに大きくなる。最小値は0、最大値は10のようである。
最小値は定義上0以上のはずだがなぜか-0.3とかが混じっていた。謎である。
比をとれば1に近い感じになった。実際の塩基配列はGCリピートなどで偏っているだろうし、変異する塩基ポジションもホットスポットのような領域があるはずなので、\frac{Ka}{Ks}が1に近いことはあまりないのではないだろうかと思う。