双序列比对的基础(2)之替换(计分)矩阵系列
主要以BLOSUM矩阵与PAM矩阵的介绍为主。 声明:该部分书中内容介绍有点少,所以我上网搜索到几篇文献和和国外大学的相关课件(从一个研究生博主处获得)。
那本篇文章就先介绍BLOSUM矩阵吧
BLOck SUBstitution Matrix:BLOSUM矩阵。 详细的来说,它们来自一组蛋白质家族中联配上的无空位区域,这些蛋白家族源于BLOCKS数据库。1
BLOSUM62矩阵广泛应用于双序列比对,也是BLAST程序默认调用的计分矩阵。
62是什么意思嘞?
BLOSUM62 is derived from Blocks containing >62%
identity in ungapped sequence alignment.2
也就是说BLOSUM62矩阵来自于序列间等同残基比超过62%的blocks(区块)。而block就是序列间联配上的无空位区域。
小小的总结 :1.先定一个阈值L,比如你最后想得到BLOSUM62矩阵,就把L定为62。2 . 前往蛋白质序列数据库,将符合序列间等同残基比大于L的序列归为一类。3.将得到的类里面的序列作多序列比对(用PAM矩阵进行的多序列对比)。4. 对比后,将保守无空位的区域划分为block。5.在block内统计频率,一个block相当于一个匹配模型。运用对数几率比得出s(a,b)。
例如:
上图的深色部分就是一个block。下面我们就以该block为例子计算残基对之间的联配的分值s(a,b)。
计算的核心就是上一篇文章提到的对数几率比(log odds ratio)即
s
(
a
,
b
)
=
log
(
p
a
b
q
a
q
b
)
s ( a , b ) = \log \left( \frac { p _ { a b } } { q _ { a } q _ { b } } \right)
s(a,b)=log(qaqbpab)。用统计得的归一化频率来代表概率。
计算流程:
c i i ( k ) c _ { i i } ^ { ( k ) } cii(k)= C n i 2 C _ { n _ { i } } ^ { 2 } Cni2
c i j ( k ) c _ { i j } ^ { ( k ) } cij(k)= n i n _ { i } ni* n j n _ { j } nj
c i j = ∑ k c i j ( k ) c _ { i j } = \sum _ { k }c _ { i j } ^ { ( k ) } cij=∑kcij(k)
c
i
j
(
k
)
c _ { i j } ^ { ( k ) }
cij(k):第k列内残基对(i,j)被观测到的次数。
n
i
n _ { i }
ni:该列中残基i被观测到的次数。
T=W* C N 2 C _ { N } ^ { 2 } CN2
q i j = c i j T q _ { i j } = \frac { c _ { i j } } { T } qij=Tcij
W:列数,N:行数。归一化的频率表示概率。
p i = q i i + ∑ j = i q i j 2 p _ { i } = q _ { i i } + \sum _ { j = i } \frac { q _ { i j } } { 2 } pi=qii+∑j=i2qij
p i p _ { i } pi:残基i在该block中出现的概率。
e i i = p i 2 e _ { i i } = p _ { i } ^ { 2 } eii=pi2
e i j = 2 p i p j e _ { i j } = 2 p _ { i } p _ { j } eij=2pipj
e i j e _ { i j } eij:残基对(i,j)随机出现的概率
s ( i , j ) = log 2 q i j e i j s(i,j)= \log_ { 2 } \frac { q _ { i j } } { e _ { i j } } s(i,j)=log2eijqij、
最后BLOSUM矩阵[i,j]=2*s(i,j),并取最邻近的整数。
本篇总结:
看到此处相信你已经对计分矩阵的建立有了一定的了解。但是,如果仔细想想会有一个疑问。直接用观测到的联配的归一化频率表示匹配模型M中的概率参数(在自然界中各残基之间的联配概率。),也就是说用一个样本的观测数据直接估计为总体的概率参数。是不是得进行最大似然估计才能估计总体的参数?是不是没进行最大似然估计?
也许看下几篇文章会得到解答。
This is a BLOSUM62矩阵.