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

MikuHatsune2013-04-15

塩基配列Complexityとしてエントロピーという統計量を使うのだが、このエントロピーの定義が自分の周辺で物議を醸している。
一般的なシャノンエントロピー
S_1=-\sum_{i\in \{A,G,C,T\}}p_i log_2 p_i
だが、前の論文では、HBVのquasispeciesとして、クローン数が全部でN個の配列で、配列s_ip_iの割合で存在しているとすると
S_2=-\sum_i^n \frac{p_i log{p_i}}{log N}
となる、と書いてある。
 
S_1は、各塩基位置において、AGCTがどの頻度で出現しているかを計算している。
S_2は、各(臨床)サンプルにおいて、ある配列をもつクローンがどの頻度で出現しているかを計算している。
と思う。
で、これらが一致するかという質問を受けたのだが、結論から言うと一致しない
S_1は、取りうる要素数C(=4)(塩基)のとき、0\leq S_1\leq log_2 Cとなる。
S_2は、未検証S_10\leq S_1\leq log_2 Cとなるので、log N(エントロピーを自然対数で計算しなおせば)で割れば0\leq S_2\leq 1になる。

library(MCMCpack)
p1 <- rdirichlet(10000, rep(1, 4))
plot(sort(rowSums(-log(p1, 2) * p1)), ylab="Entropy")


 
例えば、10塩基長、15クローンのデータがあったとして、このデータのシャノンエントロピー(S_2=0.10)となる。
各塩基でのエントロピーS_1は、総和が2.65で平均が0.265なので、どうすれば一致するんだろうか。たぶんしない。

      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
 [1,] "T"  "A"  "T"  "C"  "G"  "C"  "T"  "A"  "A"  "G"  
 [2,] "T"  "A"  "T"  "C"  "G"  "C"  "T"  "A"  "A"  "G"  
 [3,] "T"  "A"  "T"  "C"  "G"  "C"  "T"  "A"  "C"  "G"  
 [4,] "T"  "A"  "T"  "C"  "G"  "C"  "T"  "A"  "A"  "G"  
 [5,] "T"  "A"  "T"  "C"  "G"  "C"  "T"  "T"  "A"  "G"  
 [6,] "T"  "A"  "T"  "C"  "G"  "C"  "T"  "A"  "G"  "G"  
 [7,] "T"  "A"  "T"  "C"  "G"  "C"  "T"  "A"  "A"  "G"  
 [8,] "C"  "A"  "T"  "C"  "G"  "C"  "T"  "A"  "A"  "G"  
 [9,] "T"  "A"  "T"  "C"  "G"  "C"  "T"  "A"  "A"  "G"  
[10,] "T"  "A"  "T"  "C"  "G"  "C"  "T"  "A"  "A"  "G"  
[11,] "T"  "A"  "G"  "C"  "G"  "C"  "T"  "A"  "A"  "G"  
[12,] "T"  "A"  "T"  "C"  "G"  "C"  "T"  "A"  "A"  "G"  
[13,] "T"  "A"  "T"  "C"  "G"  "C"  "T"  "A"  "A"  "G"  
[14,] "T"  "A"  "T"  "C"  "G"  "T"  "T"  "A"  "T"  "G"  
[15,] "T"  "A"  "T"  "C"  "G"  "C"  "T"  "A"  "C"  "G"

library(seqinr)
library(bio3d)
bp <- 10
clone <- 15
niter <- 1000
base <- c("A", "G", "C", "T")
par(mfrow=c(1, 2))
for(m0 in c(5, 50)){
	mut_bases <- rpois(1, m0)
	res <- matrix(0, niter, 2)
	for(i0 in seq(niter)){
		s1 <- sample(base, size=bp, replace=TRUE)
		m1 <- mapply(rep, s1, clone)
		mut <- sample(seq(length(m1)), mut_bases)
		mutbase <- sample(base, size=length(mut), replace=TRUE)
		m2 <- matrix(replace(c(m1), mut, mutbase), nr=clone)
		
		#AGCTベースのエントロピー
		Sn_p1 <- mapply(function(x) table(factor(m2[ ,x], levels = base)), seq(bp)) / clone
		Sn1 <- -colSums(Sn_p1 * log(Sn_p1, 2), na.rm=TRUE)
		names(Sn1) <- seq(bp)
		#barplot(Sn1, ylim=c(0, log(length(base), 2)), xlab="base", main="Shannon entropy")
		
		#クローンベースのエントロピー
		seq0 <- table(apply(m2, 1, paste, sep="", collapse=""))
		Sn2 <- -(seq0 / clone) * log(seq0 / clone) / log(clone)
		
		res[i0, ] <- c(mean(Sn1), sum(Sn2))
	}
	xl <- yl <- range(0, res)
	plot(res, xlim=xl, ylim=yl, cex=0.7, xlab="Base", ylab="Clone")
	title(paste("Correration of Shannnon entropy\nMutation in", round(m0/(bp*clone) * 100, 1), "% in sequence"))
	abline(lm(res[, 2] ~ res[, 1]))
}

あるひとつの配列中で、変異する塩基が増えると、クローンのシャノンエントロピー(S_2)はある定数に近づくようである。

 
bio3dパッケージにあるentropy関数は、塩基エントロピー(S_1)で計算しているようである。
このなかのサンプルデータとして、HIVのプロテアーゼのどっかの領域にひっかかる。
protein BLASTより、データの1行目

PQITLWQRPLVTIKIGGQLKEALLDTGADDTVLEEMSLPGRWKPKMIGGI--GGFIKVRQYDQILIEICGHKAIGTVLVGPTPVNIIGRNLLTQIGCTLNF

を検索すると出てくる。

これのサンプルコードを実行すると、いろいろ計算して描いてくれるが、各アミノ酸残基でのエントロピーを計算している。

# Read HIV protease alignment 
aln <- read.fasta(system.file("examples/hivp_xray.fa",package="bio3d"))
# Entropy and consensus
h   <- entropy(aln)
con <- consensus(aln)
names(h$H)=con$seq
print(h$H)
H <- h$H.norm

col <- mono.colors(100)
aa  <- rev(rownames(h$freq))
oldpar <- par(no.readonly=TRUE)
layout(matrix(c(1, 2), 2, 1, byrow = TRUE), widths = 7, heights = c(2, 8), respect = FALSE)

# Plot 1: entropy
par(mar = c(0, 4, 2, 2))
barplot(H, border="white", ylab = "Entropy", space=0, xlim=c(3.7, 97.3),yaxt="n" )
axis(side=2, at=c(0.2, 0.4, 0.6, 0.8))
axis(side=3, at=(seq(0, length(con$seq), by=5) - 0.5), labels=seq(0, length(con$seq), by=5))
box()

# Plot2: residue frequencies 
par(mar = c(5, 4, 0, 2))
image(x=1:ncol(con$freq),
   y=1:nrow(con$freq),
   z=as.matrix(rev(as.data.frame(t(con$freq)))),
   col=col, yaxt="n", xaxt="n",
   xlab="Alignment Position", ylab="Residue Type")
axis(side=1, at=seq(0,length(con$seq), by=5))
axis(side=2, at=c(1:22), labels=aa, cex.axis=0.7)
axis(side=3, at=c(1:length(con$seq)), labels =con$seq, cex.axis=0.4, tick=0)
axis(side=4, at=c(1:22), labels=aa, cex.axis=0.7)
grid(length(con$seq), length(aa))
box()
for(i in 1:length(con$seq)) {
	text(i, which(aa==con$seq[i]),con$seq[i],col="white")
}
abline(h=c(3.5, 4.5, 5.5, 3.5, 7.5, 9.5, 12.5, 14.5, 16.5, 19.5), col="gray")
par(oldpar)