双序列比对的基础（2）之替换（计分）矩阵系列

  主要以BLOSUM矩阵与PAM矩阵的介绍为主。 声明：该部分书中内容介绍有点少，所以我上网搜索到几篇文献和和国外大学的相关课件（从一个研究生博主处获得）。
那本篇文章就先介绍BLOSUM矩阵吧
BLOck SUBstitution Matrix：BLOSUM矩阵。 详细的来说，它们来自一组蛋白质家族中联配上的无空位区域，这些蛋白家族源于BLOCKS数据库。¹

BLOSUM62矩阵广泛应用于双序列比对，也是BLAST程序默认调用的计分矩阵。
62是什么意思嘞？
BLOSUM62 is derived from Blocks containing >62%
identity in ungapped sequence alignment.²

  也就是说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 \left( \frac { p _ { a b } } { q _ { a } q _ { b } } \right)$。用统计得的归一化频率来代表概率。
计算流程：

$c _ { i i } ^ { ( k ) }$=$C _ { n _ { i } } ^ { 2 }$
$c _ { i j } ^ { ( k ) }$=$n _ { i }$*$n _ { j }$
$c _ { i j } = \sum _ { k }c _ { i j } ^ { ( k ) }$

 $c _ { i j } ^ { ( k ) }$:第k列内残基对（i，j)被观测到的次数。
$n _ { i }$:该列中残基i被观测到的次数。

T=W*$C _ { N } ^ { 2 }$
$q _ { i j } = \frac { c _ { i j } } { T }$

  W:列数，N：行数。归一化的频率表示概率。

$p _ { i } = q _ { i i } + \sum _ { j = i } \frac { q _ { i j } } { 2 }$

$p _ { i }$:残基i在该block中出现的概率。

$e _ { i i } = p _ { i } ^ { 2 }$
$e _ { i j } = 2 p _ { i } p _ { j }$

$e _ { i j }$:残基对（i，j）随机出现的概率

$s(i,j)= \log_ { 2 } \frac { q _ { i j } } { e _ { i j } }$、

最后BLOSUM矩阵[i,j]=2*s(i,j),并取最邻近的整数。
This is a BLOSUM62矩阵.

---

1. 《生物序列分析》 ↩︎

2. 课件地址：http://www.cs.columbia.edu/4761/assignments/assignment1/reference1.pdf ↩︎