2014-01-18

2. 配列データを圧縮する(GReEnで)：ぱらぱらめくる『Deep Sequencing Data Analysis』

大量のシークエンスデータが出るので、保管するにも移動するにも配列の圧縮は必須
その圧縮ツールの一つを紹介
算術符号(Arithmetic coding)に基礎を置く手法GReEn(入手はこちら)にある圧縮手法の一覧から
算術符号という手法
- "0110"という0/1でできた配列があったときに、それを頭から順に処理して行って、配列に特異的な線分の両端の値を出すことにする
- そうすれば、定めた線分の上の何らかの値と配列の長さを与えることで、その配列の情報を失うことなく保持できる
- 線分はどんどん短くなり、短くなればなるほど、その線分の両端を指定するための値の小数点以下桁数は増える。この桁がある意味の「情報の詳しさ」になる
- この値の桁数をなるべく増やさないためには、0/1が出たときに、出やすい方が出たらあまり短くせず、珍しい方が出たら思い切って短くするとよい
- その短くする割合をp,(1-p)で0,1に与える
- 0が配列全体でpの割合を占めるとき(と予想されるとき)、0が出たら1-p倍の長さに縮め、1が出たらp倍に縮める
- 文で説明するよりRでやる方がわかりやすい

# 算術符号処理の関数。stは0/1のベクトル。pは1が出たときに縮める割合

ar.cd.01 <- function(st,p){
	ret <- c(0,1)
	ret.hx <- matrix(0,length(st)+1,2)
	ret.hx[1,] <- ret
	for(i in 1:length(st)){
		if(st[i]==0){
			ret[2] <- ret[2] - (ret[2]-ret[1])*(1-p)
		}else{
			ret[1] <- ret[1] + (ret[2]-ret[1])*p
		}
		ret.hx[i+1,] <- ret
	}
	return(list(ans=ret,ans.hx=ret.hx))
}

# ある程度のルールのある0/1列を作る
st <- c()

n.st <- 20
st <- c(0)
for(i in 2:n.st){
	if(st[i-1]==0){
		if(runif(1) > 0.3){
			st[i] <- 0
		}else{
			st[i] <- 1
		}
	}else{
		if(runif(1) > 0.4){
			st[i] <- 1
		}else{
			st[i] <- 0
		}
	}
}
# できた配列を見る
st
# 1の割合を見る
p <- sum(st)/length(st)
# 算術符号処理
out <- ar.cd.01(st,p)
# 線分がどのように縮んでいくかを図示
matplot(out$ans.hx,type="l")

[f:id:ryamada22:20140117135649j:image]
> st
 [1] 0 0 0 0 0 1 1 0 0 0 1 1 1 1 1 0 0 0 1 1

- これを0/1ではなくて、もっとたくさんの文字にすると…
  - $X=(x_1,x_2,...,x_k)$ なるk種類の文字列があって、それぞれが出たときの「縮み率」を $P=(p_1,...,p_k);\sum_{i=1}^k p_i=1$ とすると
  - この $P$ を使って $Q=(q_1=0,q_2=p_1,p_1+p_2,p_1+p_2+p_3,...,q_{u+1}=\sum_{i=1}^u p_i,...,q_{k+1}=\sum_{i=1$ ^k p_i=1)]なる長さ $k+1$ のベクトルが作れて、
  - $[a,b);0\le a < b \le 1$ なる線分に対して、 $x_i$ が出たら $[a*(1-q_i)+b*q_i,a*(1-q_{i+1})+b*q_{i+1})$ に変更するという処理をすればよい
  - これもRでやっておこう

ar.cd <- function(st,p){
	n <- length(p)
	q <- c(0,cumsum(p))
	ret <- c(0,1)
	ret.hx <- matrix(0,length(st)+1,2)
	ret.hx[1,] <- ret
	for(i in 1:length(st)){
		tmp <- st[i] + 1
		tmp.p <- q[tmp:(tmp+1)]
		tmp.p.2 <- 1-tmp.p
		ret[1] <- ret[1]*tmp.p.2[1]+ret[2]*tmp.p[1]
		ret[2] <- ret[1]*tmp.p.2[2]+ret[2]*tmp.p[2]
		ret.hx[i+1,] <- ret
	}
	return(list(ans=ret,ans.hx=ret.hx))
}

st <- sample(0:3,4,replace=TRUE)
st
p <- runif(4)
p <- p/sum(p)

out2 <- ar.cd(st,p)
out2
# 線分がどのように縮んでいくかを図示
matplot(out2$ans.hx,type="l")

> out2
$ans
[1] 0.9786711 0.9809776

$ans.hx
          [,1]      [,2]
[1,] 0.0000000 1.0000000
[2,] 0.5482826 0.9885911
[3,] 0.7896960 0.9835676
[4,] 0.9786711 0.9835676
[5,] 0.9786711 0.9809776

さて、GReEnはどうやっているかというと：
- ATGCの４文字があるので、0/1ではなくて4文字バージョン
- どこが工夫か、というと、「どの文字が出たら」どのくらいの縮み率にするかのところ
- Reference配列があれば、「referenceと同じことが多い」はずなので、そのときは縮み率が小さいようにして、それ以外のときに縮むようにすると、圧縮効率がよい。これを用いて、塩基箇所ごとに圧縮率を調整する、という工夫を入れている(copy model)
- また、Reference配列に「合っている」とみなして縮み率を決めるとは言うものの、アラインメントして「refecence配列の何番目の塩基」と「処理しているtarget配列の何番目の塩基」とが対応関係にあるかは、ずれて行くので、基準を定めて、その「合っている位置」の位置取りのやり直しをする。その上で、同じ位置取りの状況で１塩基ずつ処理していって、「良く合った(readで言えば、大文字のATGC(精度のよい塩基コール)として合致」と「一応、合った(readで言えば、小文字のatgc(精度がいまいちの塩基コール)として合致」との情報をとり、それを「縮み率」に反映させながら進める(合致がよいところは、referenceと同じ塩基が出る確率が高いが、合致率が悪ければ、そうではない、という情報を使う)
- ｢位置の取り直し」をした後、n塩基の処理が終わったとする。そこまでの処理でa塩基は塩基も正しく、大文字だった、b塩基は塩基は正しかったが小文字だったとする。もちろん $n \ge a+b$
- このときに、次の塩基にどういうPを想定しているか、というと、 $\frac{a+1}{n+3},\frac{b+1}{n+3},\frac{n-(a+b)+1}{n+3}$ の割合で、「大文字の当たり(referenceの塩基の大文字)」「小文字の当たり」「外れ」とする
- この分子・分母に出てくる1,3というのは、ディリクレ分布で出てくる「補正」の値。たとえば、a=b=0のように情報がない場合には、大文字がよいか小文字がよいか半々にするのがよいし…という話。 $(x_1,x_2,...,x_k)$ の観察があったときに、 $\frac{x_j+1}{\sum_{i=1}^k x_i + k}$ が期待値であるということ(これは二項情報からのベータ分布推定、多項情報からのディリクレ分布推定のことであるし、この本では、"Laplace–Bayes estimator"として説明してある
- ただし、「外れ」には、referenceとは違う複数の塩基があって、その大文字と小文字とがあるので、その割合も決めないといけない。それは、どの塩基の大文字/小文字がtarget配列に占める割合で重みづけ分配している