BLOSUM62是应用得最广的氨基酸替换矩阵，是BLAST中蛋白质序列比对的默认矩阵。
1. 如何下载替换计分矩阵
ftp://ftp.ncbi.nih.gov/blast/matrices/
BLOSUM62
https://www.ncbi.nlm.nih.gov/Class/FieldGuide/BLOSUM62.txt
ftp://ftp.ncbi.nih.gov/blast/matrices/BLOSUM62
2. 推导过程
2.1 blocks的发展
BLOSUM替换计分矩阵由Steven Henikoff和Jorja Henikoff 在1992年提出，是一个log-odds 的矩阵，这里的odds指什么（下文有）。
在蛋白质序列数据库中，基于序列之间的identity（大于一个阈值）将这些蛋白质序列cluster为500个group，每个group里面的序列做多序列比对，将保守无gap的区域划分为block，一共2000多个。identity可以取很多值，BLOSUM62矩阵用的是identity大于等于62%。然后基于这些block，可以找出20种aa所有替换情况（210）的频率。
一个block如下：
到1996年，基于770蛋白质家族，已有大约3000个blocks
如何理解：Substitution frequencies weigh more heavily by protein sequences having less than 62% identity. Therefore, BLOSUM62 is useful for aligning and scoring proteins that show less than 62% identity.
2.2 计算过程
the (log) ratio of the observed probability of substitution of one amino acid by another divided by the probability expected purely due to chance。
翻译过来就是：一个氨基酸被另一个氨基酸替换所观察到的频率除以它俩因为随机而出现在一起的概率，然后取log值。
随机模型R：每一个碱基都是以频率q独立出现
匹配模型M：匹配上的碱基对以联合概率p(ab)出现，p(ab)的值可以认为是b就是由a演变而来的概率
定义
2.2.1 计算出一个block中每一列每一种配对出现的频数
以上文的block为例，第一列为AABACA，分别求出可能出现的配对情况（两两配对）和次数：
| 配对 | 次数 | 公式 |
|---|---|---|
| AA | 6次 | 4*3/2 |
| AB | 4次 | 4*1 |
| AC | 4次 | 4*1 |
| BB | 0次 | 1*0/2 |
| BC | 1次 | 1*1 |
| CC | 0次 | 1*0/2 |
于是可以归纳出一般计算配对频数的方法：若碱基相同，n*(n-1)/2；若碱基不同，n1*n2，n为碱基出现次数。
2.2.2 遍历block的每一列，将特定配对情况的频数都加起来
以AB配对为例，在第二列中AB的频数为8次，其他列没有了，所以一共是12次。同理：
| 配对 | 频数 |
|---|---|
| AA | 7 |
| AB | 12 |
| AC | 9 |
| AD | 5 |
| BB | 26 |
| BC | 11 |
| BD | 0 |
| CC | 20 |
| CD | 5 |
| DD | 10 |
这些频数的加和一定是等于这个block中所有能观察到的配对数的和：
其中，w为列数，n为行数
2.2.3 计算出每一种配对情况的观察到的频率
以AB配对为例：
2.2.4 基于block计算某种氨基酸出现的概率
以A为例，AA配对贡献两个A，A(其他氨基酸)这类配对贡献一个，氨基酸在配对过程中总的出现次数是105*2
拓展到其它情况，氨基酸i出现的概率：
2.2.5 计算由于随机因素两个氨基酸一起出现的期望概率
2.2.6 求出log odds ratio
2.2.7 确定计分矩阵
上面的值再乘以2四舍五入取整即可
3. PAM and BLOSUM的异同点
- PAM矩阵的构建基于进化模型——需要一个突变率，其估计需要通过构建系统发生树及推断祖先序列。BLOSUM矩阵的构建如上面讲的，替换率可以直接从保守无gap的block中看出来。因此PAM经常被用来构建系统发生树，BLOSUM经常被用来做局部比对。
- PAM矩阵的构建是对全长序列（包括保守序列以及不保守序列）进行全局比对，而BLOSUM矩阵对block进行局部比对。通过比较BLOSUM62和PAM160可以发现，BLOSUM62对于亲水氨基酸的替换更严格，对于疏水氨基酸的替换更容忍。此外，对于罕见的氨基酸，比如半胱氨酸和色氨酸，BLOSUM也更容忍。
- 适用范围
PAM矩阵
后缀数字越大，处理进化上距离远的蛋白质序列效果越好；
后缀数字越小，处理进化上距离近的蛋白质序列效果越好。
BLOSAM矩阵
后缀数字越大，处理进化上距离近的蛋白质序列效果越好。
如果事先并不知道序列之间的相似程度，可以多试几个矩阵，最后再使用给分最高的那个矩阵
4. 参考
两篇古老的文献
http://sci-hub.se/10.1073/pnas.89.22.10915
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC329220/pdf/nar00103-0200.pdf
一点章节内容
https://www.sciencedirect.com/topics/medicine-and-dentistry/blosum
一份PPT
http://www.cs.columbia.edu/4761/assignments/assignment1/reference1.pdf
一本书
《生物序列分析》王俊等人翻译
网友评论