【5.4】Position-Specific Iterated (PSI)-BLAST

PSI-BLAST(Position-Specific Iterated BLAST)是很灵敏的BLAST工具，可以用来找距离很远的蛋白或者找蛋白家族新的成员。

一、psi-blast的意义

图解PSI-Blast的意义

注： 这个PSMM怎么用，后续还得进一步了解一下。。。

二、psi-blast原理与操作示范

第一步的psiblast其实就是简单的blast,构建新的矩阵，然后在对于大于阈值重新blast,比对的结果中有用黄色标记的就是之前没有找序列。

为了提高速度，标准 BLAST 牺牲了一定的准确度，牺牲掉的准确度对高度相似的序列， 也就是亲缘关系近的序列构成不了威胁，不会把它们落掉，但是对于那些只有一点点相似， 也就是远源的序列，就有点麻烦了。它们很可能被落掉而没有被 BLAST 发现。

要解决这个问题，可以用 PSI-BLAST。PSI-BLAST 可以搜罗出一个庞大的蛋白质家族， 当然也包括标准 BLAST 不小心漏掉的那些远房亲戚。换言之，标准 BLAST 找到了直接认识 的朋友，但朋友的朋友都丢掉了。PSI 是 Position-Specific Iterated 首字母缩写，中文是位点特 异性迭代。PSI-BLAST 的特色是搜完一遍再搜一遍，且从第二次搜索开始，每次搜索前都利 用上一次搜索到的结果创建一个位置特异权重矩阵以扩大本次搜索的范围。如此反复直至没 有新的结果产生为止。位置特异权重矩阵(Position-Specific Scoring Matrix，简称 PSSM)是 以矩阵的形式，统计一个多序列比对中，每个位置上不同残基出现的百分比。假设 A 的朋 友只有 B，B 朋友除 A 外还有 C。如果输入序列的第一个位置是 A，那么在第一轮没有 PSSM 辅助的情况下，只有第一个位置是 A 或 B 的序列被找到了。它们是图 1-A 中所示的四条序 列。根据这四条序列创建的 PSSM(图 1-B)得知，第一个位置可以是 A，也可以是 B，那么 在第二轮搜索中，除了 A 的朋友 B 之外，B 的朋友 C 也可以出现在第一个位置了。这样如 此反复，我们就可以找到朋友的朋友了。

第一轮搜索结果(图 2)和标准 BLAST 是一样的。略有不同的是，在第三部分(Descriptions) 列表的最右侧多了两列。绿色标记的一列是将要为下一轮搜索创建 PSSM 所选择的序列，默 认第一轮搜索到的序列都将用来创建 PSSM。红色标记的一列是已在本轮搜索中用来创建 PSSM 的序列。因为是第一轮搜索，之前还没有搜索到任何序列，也就是第一轮的搜索过程中没有使用 PSSM，所以这一列都为空。接下来点“go”进行第二轮搜索。“go”左侧的输入 框里可以指定列出搜索到的前多少条序列。因为 PSI-BLAST 无休止的一直做下去，也许会 把数据库中的全部序列都找出来。

第二轮的搜索结果里(图 3)，红色标记的一列，也就是已在第二轮搜索中被用来创建 PSSM 的序列，大部分都打勾了。但是，再往后面会看到有些标黄的序列没有打勾。这些没 有打勾的标黄的序列就是在第二轮搜索中新找到的序列。它们将用于创建下一轮搜索使用的 PSSM，但是在本轮搜索中，它们没有被用到，所以没有打勾。而没有标黄的这些打勾的序 列，是在第一轮中就已经找到的序列。

仔细看一下这些新找到的序列，它们和输入序列之间的一致度其实并不比第一轮中找到 的序列低。但是因为算法上的原因，它们被标准 BLAST 漏掉了。好在通过 PSI-BLAST，它 们又被找回来了。

如果想要快速找到本轮搜索新找到的序列，也就是第一条标黄的序列，可以点击结果页 面上方“skip to the first new sequence”链接(图 5)。

三、PSSM详解

示例：

cmd_2 = 'psiblast -query %s -evalue 0.0001 -db %s -num_iterations 3 -out_ascii_pssm %s  -out  %s -out_pssm  %s  -num_threads %s' % \
        (input_fa, blast_result_db, pssm_result, pssm_blast_result, out_pssm_result,num_threads) 

参数解释：

• query 输入的fasta序列
• db 比对的数据库
• num_iterations 用构建的打分矩阵，重新搜索的次数（次数越多，越能找到离得远的物种）
• out_ascii_pssm pssm打分矩阵，前半部分是打分矩阵，后半部分是每个位置对应的残基比例
• out 详细的比对结构，json数据格式
• out_pssm 比对序列位置的概况
• num_threads 线程数

#获取比对结果中，相似度大于30%的序列

identity_cutoff = 30
cmd = """  awk 'BEGIN{tmp[000]=1} {if($3>%s && (!($2 in tmp)) ){print $0;tmp[$2]=1}}' %s > %s """ % (identity_cutoff,blast_result_raw,blast_result)
os.system(cmd)

参考资料

• 山东大学 生物信息学课题组荣誉出品 巩晶老师课件
• https://www.pharm.auth.gr/papadopoulou/students/bioinformatics/2006-09-13_lect04_ch5.pdf