バイオインフォマティクス

生命情報学に関する基本的な理論や解析ソフトウェアの使用方法について.

ダイナミックプログラミング (DP) は,そのままの形式で解くには複雑すぎる問題を小さい簡単な問題に分割して解を導き出すための方法.動的計画法とも呼ばれる.バイオインフォマティクスの分野では,配列を全長に渡り無理やり揃える方法であるグローバルアライメント法 (Needleman-Wunsch法),配列間のとてもよく似ている部分だけのアライメントを行うローカルアライメント法 (Smith-Waterman法),グローバルアライメントにおいて開始と末端におけるギャップにペナルティを課さないグローカル法において活用されている.

グローバルアライメントの実演

ここでは,MICKEY と MINNY というふたつのアミノ酸配列を揃えたいとする.このとき,最初に以下のようなマトリックスを生成する.これをDP行列 ($H$) と呼ぶ.6文字のMICKEYに対して7行,5文字のMINNYに対して6列の行列となる.以降,一般化のために最初の配列長さをM,もう1本の配列の長さをNとする.

MINNY
012345
0
M1
I2
C3
K4
E5
Y6

最初に,ギャップペナルティの値を0列および0行のセルに対して計算する.生物の進化の過程において,元々相同な関係にあった配列間でアミノ酸に置換や欠損が生じることがあるが,それらの座位の片方にはその置換か欠損をを意味するギャップという文字 - を入れる.それがアライメントの結果に入ることに対するペナルティのことをギャップペナルティという.アミノ酸からアミノ酸への変異より欠損とか挿入等のイベントはより起こりにくいだろうという仮定のもと,ペナルティが課されている.オープンギャップペナルティ ($\mu$) を-10,拡張ギャップペナルティ ($\sigma$) を-2とすると以下のようにDP行列は埋められる.

MINNY
012345
00-10-12-14-16-18
M1-10
I2-12
C3-14
K4-16
E5-18
Y6-20

次に,各セルを埋めていく.このとき行 ($i$) のアミノ酸と列 ($j$) アミノ酸が揃う時にそのアミノ酸に与えるスコアを利用する.以下のようなスコアテーブルである.これをアミノ酸置換行列と呼ぶ.他にもたくさんの種類があるが,ここで用いるのは最も良く用いられている行列,BLOSUM62 である.あるアミノ酸が別のアミノ酸へどれくらいの変異し易さを持つかということを示す行列であり,アミノ酸とアミノ酸の交差点の値が小さいほど,それらのアミノ酸間で変異が起き易いことを表す.

   A  R  N  D  C  Q  E  G  H  I  L  K  M  F  P  S  T  W  Y  V  B  Z  X  *
A  6 -2 -2 -3 -1 -1 -1  0 -2 -2 -2 -1 -1 -3 -1  2  0 -4 -3  0 -2 -1 -1 -6
R -2  8 -1 -2 -5  1  0 -3  0 -4 -3  3 -2 -4 -3 -1 -2 -4 -3 -4 -2  0 -2 -6
N -2 -1  8  2 -4  0  0 -1  1 -5 -5  0 -3 -4 -3  1  0 -6 -3 -4  5  0 -2 -6
D -3 -2  2  9 -5  0  2 -2 -2 -5 -5 -1 -5 -5 -2  0 -2 -6 -5 -5  6  1 -2 -6
C -1 -5 -4 -5 13 -4 -5 -4 -4 -2 -2 -5 -2 -4 -4 -1 -1 -3 -4 -1 -5 -5 -3 -6
Q -1  1  0  0 -4  8  3 -3  1 -4 -3  2 -1 -5 -2  0 -1 -3 -2 -3  0  5 -1 -6
E -1  0  0  2 -5  3  7 -3  0 -5 -4  1 -3 -5 -2  0 -1 -4 -3 -4  1  6 -1 -6
G  0 -3 -1 -2 -4 -3 -3  8 -3 -6 -5 -2 -4 -5 -3  0 -2 -4 -5 -5 -1 -3 -2 -6
H -2  0  1 -2 -4  1  0 -3 11 -5 -4 -1 -2 -2 -3 -1 -3 -4  3 -5 -1  0 -2 -6
I -2 -4 -5 -5 -2 -4 -5 -6 -5  6  2 -4  2  0 -4 -4 -1 -4 -2  4 -5 -5 -2 -6
L -2 -3 -5 -5 -2 -3 -4 -5 -4  2  6 -4  3  1 -4 -4 -2 -2 -2  1 -5 -4 -2 -6
K -1  3  0 -1 -5  2  1 -2 -1 -4 -4  7 -2 -5 -2  0 -1 -4 -3 -3 -1  1 -1 -6
M -1 -2 -3 -5 -2 -1 -3 -4 -2  2  3 -2  8  0 -4 -2 -1 -2 -1  1 -4 -2 -1 -6
F -3 -4 -4 -5 -4 -5 -5 -5 -2  0  1 -5  0  9 -5 -4 -3  1  4 -1 -5 -5 -2 -6
P -1 -3 -3 -2 -4 -2 -2 -3 -3 -4 -4 -2 -4 -5 11 -1 -2 -5 -4 -4 -3 -2 -2 -6
S  2 -1  1  0 -1  0  0  0 -1 -4 -4  0 -2 -4 -1  6  2 -4 -3 -2  0  0 -1 -6
T  0 -2  0 -2 -1 -1 -1 -2 -3 -1 -2 -1 -1 -3 -2  2  7 -4 -2  0 -1 -1 -1 -6
W -4 -4 -6 -6 -3 -3 -4 -4 -4 -4 -2 -4 -2  1 -5 -4 -4 16  3 -4 -6 -4 -3 -6
Y -3 -3 -3 -5 -4 -2 -3 -5  3 -2 -2 -3 -1  4 -4 -3 -2  3 10 -2 -4 -3 -2 -6
V  0 -4 -4 -5 -1 -3 -4 -5 -5  4  1 -3  1 -1 -4 -2  0 -4 -2  6 -5 -4 -1 -6
B -2 -2  5  6 -5  0  1 -1 -1 -5 -5 -1 -4 -5 -3  0 -1 -6 -4 -5  5  0 -2 -6
Z -1  0  0  1 -5  5  6 -3  0 -5 -4  1 -2 -5 -2  0 -1 -4 -3 -4  0  5 -1 -6
X -1 -2 -2 -2 -3 -1 -1 -2 -2 -2 -2 -1 -1 -2 -2 -1 -1 -3 -2 -1 -2 -1 -2 -6
* -6 -6 -6 -6 -6 -6 -6 -6 -6 -6 -6 -6 -6 -6 -6 -6 -6 -6 -6 -6 -6 -6 -6  1

この置換行列を用いて,以下のような計算式を考る.これらから,$d$,$u$ および $l$ を計算する.ここで,$H$ はDP行列を示す.また,$I$ と $J$ も行列で,これはどちらもM行N列,すなわち揃えるアミノ酸の長さに由来する要素からなる行列.このふたつの行列はあらかじめものすごく小さい値 (-999999とか) で初期化されているものとする.また,$a_i$ と $b_j$ はそれぞれ,1本目と2本目の配列の $i$ 番目および $j$ 番目のアミノ酸.$s()$ はアミノ酸置換行列を表す.

\begin{eqnarray*}d=\max \begin{Bmatrix} H_{i-1,j-1}+s(a_i,b_j) \\ I_{i-1,j-1}+s(a_i,b_j) \\ J_{i-1,j-1}+s(a_i,b_j) \\ \end{Bmatrix}\tag{1}\end{eqnarray*}
\begin{eqnarray*}u=\max \begin{Bmatrix} H_{i-1,j}+\mu \\ I_{i-1,j}+\sigma \\ \end{Bmatrix}\tag{2}\end{eqnarray*}
\begin{eqnarray*}l=\max \begin{Bmatrix} H_{i,j-1}+\mu \\ J_{i,j-1}+\sigma \\ \end{Bmatrix}\tag{3}\end{eqnarray*}

次に,以下の式から $H$ の値を計算する.この計算を繰り返してDP行列 $H$ の全てのセルを埋める.またこのとき,上の計算式において $H_{i,j}$ が $(d,u,l)$ のどれに由来するものかを覚えておかなければならない.

\begin{eqnarray*}H_{i,j}=\max \begin{Bmatrix} d,u,l\\ \end{Bmatrix}\tag{4}\end{eqnarray*}

これらは,新たに $P$ というポインタ行列を用意して,そこに記録する.これは,DP行列と同じく(M+1)行(N+1)列の行列となる.初期化は以下のように行う.

MINNY
012345
0-LLLLL
M1U
I2U
C3U
K4U
E5U
Y6U

実演する.ここでは,$(i,j)=(1,1)$ を求めるとすると式は以下のようになる.

\begin{eqnarray*}d=\max \begin{Bmatrix} H_{0,0}+s(M,M) &=& 0+8\\ I_{0,0}+s(M,M) &=&-999999+8 \\ J_{0,0}+s(M,M) &=&-999999+8\\ \end{Bmatrix}=8\tag{5}\end{eqnarray*}
\begin{eqnarray*}u=\max \begin{Bmatrix} H_{0,1}-10 &=&-10-10\\ I_{0,1}-2&=&-999999-2 \\ \end{Bmatrix}=-20\tag{6}\end{eqnarray*}
\begin{eqnarray*}l=\max \begin{Bmatrix} H_{1,0}-10 &=&-10-10\\ J_{1,0}-2 &=&-999999-2\\ \end{Bmatrix}=-20\tag{7}\end{eqnarray*}

これより,$H_{1,1}$ は以下のように計算される.また,その値は $d=8$ に由来する.

\begin{eqnarray*}H_{1,1}=\max \begin{Bmatrix} 8,-20,-20\\ \end{Bmatrix}=8\tag{8}\end{eqnarray*}

よって,DP行列は以下のように更新される.

MINNY
012345
00-10-12-14-16-18
M1-108
I2-12
C3-14
K4-16
E5-18
Y6-20

また,ポインタ行列は以下のように更新される.DP行列で埋めた8という値がdに由来するものだからDと書くことになる.

MINNY
012345
0-LLLLL
M1UD
I2U
C3U
K4U
E5U
Y6U

他の全てのセルについて同じ計算をして以下のようにDP行列を完成させる.

MINNY
012345
00-10-12-14-16-18
M1-108-2-4-6-8
I2-12-214420
C3-14-44100-2
K4-16-624100
E5-18-80247
Y6-20-10-2-3-314

また,ポインタ行列は以下のように完成される.

MINNY
012345
0-LLLLL
M1UDLLLL
I2UUDLLL
C3UUUDDD
K4UUUDDL
E5UUUDDD
Y6UUUDDD

最後に,トレースバックという作業を行う.最初に,ポインタ行列の最も右かつ最も下のセル,$(i,j)=(6,5)$ の値を取得する.この場合,D となる.

D

次に,この値が D だったので左上のセルの値を取得する.この D はdiagonalを意味する D.また,U は up,L は left を意味する.もし,セルの値を取得して,その値が D なら次は左上のセルを,U なら上のセルを,L なら左のセルの値を取得する.この場合左上のセルの値は D なので以下のようになる.

DD

以上のような手順で一番左上のセルに到達するまでの値を取得する.以下のようになる.

DDDUDD

次にこの文字に沿ってアライメントを行う.D は2本の配列の文字が揃ったことを,L は1本目の配列にギャップ - が入ったことを,U は2本目の配列にギャップが入ったことを意味する.1本目の配列,MICKEY を逆順にした YEKCIM と2本目の配列を逆順にした YNNIM を揃える.最初の文字は D なので最初 (逆順にする前の最後) の文字はアラインされることになる.

Y
|
Y

次もDなので,アラインされる.

YE
||
YN

繰り返す.

YEK
|||
YNN

次の文字は U なので,この場合,2本目の配列にギャップが入る.以下のようになる.

YEKC
||| 
YNN-

同じように繰り返す.

YEKCI
||| |
YNN-I

最後の文字まで到達したのでトレースバックは終了する.

YEKCIM
||| ||
YNN-IM

最後に,アライメントを元の順に戻してアライメントは終了.これがダイナミックプログラミングとBLOSUM62というアミノ酸置換行列を用いて計算された最適解となる.

MICKEY
|| |||
MI-NNY

Pythonによる各種アライメント法の実装

アライメントのアルゴリズムには上のグローバルアライメント法以外にローカルアライメント法とセミグローバルアライメント法がある.これらも併せて実装したのが以下のプログラム.このとき,オプションのフラッグ -a を追加し,その引数に nw を指定したらグローバルアライメント法が,sw を指定したらローカルアライメント法が,sg を指定したらセミグローバルアライメント法が実行される.

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
import argparse,re

parser=argparse.ArgumentParser(description="This program aligns two sequences by smith-waterman, needleman-whunch or semi-global alignment method.")
parser.add_argument('-i',type=str,dest="input",required=True,help='Input pairwise fasta file name.')
parser.add_argument('-s',type=str,dest="matrix",required=True,help='Amino acid substitution matrix file name.')
parser.add_argument('-f',type=int,dest="op",required=True,help='Gap open penalty. Normally, negative integer.')
parser.add_argument('-g',type=int,dest="ep",required=True,help='Gap extension penalty. Normally, negative integer.')
parser.add_argument('-a',type=str,dest="mode",default="sg",help='Algorithm. [nw|sg|sw]. Default valude is \'sg\'.')
args=parser.parse_args()

def main():
	name,seq=readfasta(args.input)
	alseqi,alseqj=align(seq[0],seq[1],args.op,args.ep,args.matrix,args.mode)
	print(name[0],"\n",alseqi,"\n",name[1],"\n",alseqj,"\n",sep="",end="")
	
def align(seqi,seqj,op,ep,matrix,mode):
	score=readmatrix(matrix)
	NONE,LEFT,UP,DIAG=0,1,2,3
	maxi,maxj=len(seqi),len(seqj)
	H=[[0 for j in range(maxj+1)] for i in range(maxi+1)]
	I=[[-999999 for j in range(maxj+1)] for i in range(maxi+1)]
	J=[[-999999 for j in range(maxj+1)] for i in range(maxi+1)]
	P=[[0 for j in range(maxj+1)] for i in range(maxi+1)]
	
	if mode=="sg" or mode=="nw":
		for i in range(len(P)):
			for j in range(len(P[0])):
				if i==0 and j==0:
					pass
				elif i==0:
					P[i][j]=LEFT
				elif j==0:
					P[i][j]=UP
		if mode=="nw":
			for i in range(len(H)):
				for j in range(len(H[0])):
					if i==0 and j==0:
						pass
					elif i==0:
						H[i][j]=op+(j-1)*ep
					elif j==0:
						H[i][j]=op+(i-1)*ep
	
	pi,pj=0,0
	for i in range(1,maxi+1):
		ai=seqi[i-1]
		for j in range(1,maxj+1):
			aj=seqj[j-1]
			
			I[i][j]=max(H[i-1][j]+op,I[i-1][j]+ep)
			J[i][j]=max(H[i][j-1]+op,J[i][j-1]+ep)
			diagscore=max(H[i-1][j-1]+score[ai][aj],I[i-1][j-1]+score[ai][aj],J[i-1][j-1]+score[ai][aj])
			leftscore=J[i][j]
			upscore=I[i][j]
			maxscore=max(diagscore,upscore,leftscore)
			
			H[i][j]=maxscore
			if mode=='sw':
				H[i][j]=max(0,maxscore) 
			
			if mode=='sw':
				if H[i][j]==0:
					pass
				elif maxscore==diagscore:
					P[i][j]=DIAG
				elif maxscore==upscore:
					P[i][j]=UP
				elif maxscore==leftscore:
					P[i][j]=LEFT
			else:
				if maxscore==diagscore:
					P[i][j]=DIAG
				elif maxscore==upscore:
					P[i][j]=UP
				elif maxscore==leftscore:
					P[i][j]=LEFT
			pj=j
		pi=i
	
	if mode=="sw":
		pi,pj=findmaxindex(H)
	elif mode=="sg":
		rowmax,rowmaxindex=findmaxvalueindex(H[-1])
		tmp=[]
		for i in range(len(H)):
			for j in range(len(H[i])):
				if j==maxj:
					tmp.append(H[i][j])
		colmax,colmaxindex=findmaxvalueindex(tmp)
		if rowmax>colmax:
			for k in range(len(P[maxi])):
				if k>rowmaxindex:
					P[maxi][k]=LEFT
		else:
			for k in range(len(P)):
				if k>colmaxindex:
					P[k][maxj]=UP
	
	pv=P[pi][pj]
	alseqi,alseqj="",""
	while pv!=NONE:
		if pv==DIAG:
			pi-=1
			pj-=1
			alseqi+=seqi[pi]
			alseqj+=seqj[pj]
		elif pv==LEFT:
			pj-=1
			alseqi+='-'
			alseqj+=seqj[pj]
		elif pv==UP:
			pi-=1
			alseqi+=seqi[pi]
			alseqj+='-'
		pv=P[pi][pj]
	alseqi=alseqi[::-1]
	alseqj=alseqj[::-1]
	
	return [alseqi,alseqj]

def findmaxindex(ls):
	maxindex=[0,0]
	max=-999999
	for i in range(len(ls)):
		for j in range(len(ls[i])):
			if max<ls[i][j]:
				max=ls[i][j]
				maxindex[0],maxindex[1]=i,j
	return maxindex

def findmaxvalueindex(ls):
	maxindex=0
	max=-999999
	for i in range(len(ls)):
		if max<ls[i]:
			max,maxindex=ls[i],i
	return [max,maxindex]

def readfasta(input):
	name,seq=[],[]
	fin=open("%(input)s"%locals(),'r')
	tmpline=""
	for line in fin:
		line=line.rstrip()
		if line.startswith('>'):
			name.append(line)
			if tmpline:
				seq.append(tmpline)
				tmpline=""
		else:
			tmpline+=line.upper()
	seq.append(tmpline)
	fin.close()
	if len(name)!=2 or len(seq)!=2:
		sys.exit("Input sequence file must have just two sequences.")
	return [name,seq]
	
def readmatrix(input):
	matrix={}
	fin=open("%(input)s"%locals(),'r')
	head=0
	aa=[]
	for line in fin:
		line=line.strip()
		if not line.startswith('#'):
			if head==0:
				aa=re.split('\s+',line)
				head+=1
			else:
				tmp=re.split('\s+',line)
				tmpdict={}
				rowname=tmp.pop(0)
				for i in range(len(tmp)):
					tmpdict[aa[i]]=int(tmp[i])
				matrix[rowname]=tmpdict
	fin.close()
	return matrix

if __name__ == '__main__':
	main()
Hatena Google+