blast--remove the redundant protein chains whose sequence identity ≤30%

半学期在研究蛋白质相互作用位点的预测，由于生物短板，所以遇到很多困难。不说这些心酸事了，说一下这篇博客的内容吧。我仔细研读了一篇关于蛋白质相互作用位点预测的论文：Predicting protein interaction sites from residue spatial sequence profile and evolution rate（Bing Wanga,b,1, Peng Chena,b,1, De-Shuang Huanga,*, Jing-jing Lia,Tat-Ming Lokc, Michael R. Lyud）这篇论文的数据是来自Fariselli et al研究的113对蛋白质链，对着226条蛋白质链进行了去除冗余，使用了blast软件，本文主要就介绍如何使用此软件实现去除冗余的过程。原文如下：

After removal of redundant chains, we obtained a data set of 69 protein
chains (sequence identity <30%); all the data are available upon request.

Blast安装

在ftp://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/LATEST/ 下载对应版本，我下载的是windows版本，直接安装即可，添加环境变量即可使用。

安装完成后，在安装目录里可以看到很多功能，比如blastp,makeblastdb就是我所使用的两个程序。如图

创建自定义数据库makeblastdb

需要做的是在226条蛋白质链中，从几条很是相似的蛋白质链中只保留一条，剔除其余几条，然后就可以得到一些相似度均小于30%的蛋白质链数据集。blast比对主要是一对一，比对蛋白质序列的一致性，同源性，我们需要把一条蛋白质链与其他225条蛋白质链进行比对。首先，需要创建一个226条链的蛋白质链数据集，实现一对多比对。
blast中有makeblastdb程序，可以自主创建数据集。具体代码是：

import subprocess#这个包，可以实现执行命令行subprocess.call(['makeblastdb','-in','pplist.fasta','-dbtype','prot','-title','pplist','-out','pplist'])#-in 后面跟所有序列的综合文件#-dtype 后面跟这个数据库的类型，prot是蛋白质#-title 后面跟，你给这个数据库起的名字

你的数据库就建好了，名字是pplist,下文进行blastp时会使用到。

蛋白质序列比对blastp

Blastn 核酸-核酸比对

Blastp 蛋白-蛋白比对

Blastx 核酸-蛋白比对

使用blastp进行蛋白质比对，提取比对信息中的iden信息，生成一个226*226的矩阵。
blastp使用方法：

subprocess.call(['blastp','-db','pplist','-query' ,'E:\Users\Penny\PycharmProjects\BPNNproteinsites\\fasta\\'+ pdbchain + '.fasta ','-outfmt','6','-out','blastp1\\'+ pdbchain +'.txt'])#-db 后面跟你创建的数据库，如果比对的是非冗余或者其他数据库也可以，根据需要决定#-query 后面跟你要比对的序列（文件）#-outfmt 后面跟数据输出的格式，一共有12中，我使用的是6，表格的形式#-out 后面跟的是结果输出的文件，如果没有此项，将输出在屏幕上

这里是命令行显示的结果，目前我还不是很明白，为什么每条比对有两条信息，请高手指点!其中的第三列信息，就是我们所需要的iden信息，如果比对的结果大于30%，我们将iden这个矩阵中这个位点的值设为1，否则，设置为0.
这里，还有一个问题，就是显示的结果，并不是数据库中所有的蛋白质链都有比对结果，我认为是两条链的相关性，不大，我将其位点设置为0.
最终我生成一个226*227的矩阵，矩阵的第一列用来保存蛋白质链的信息，方便后来代码的识别。

iden矩阵进行筛选蛋白质链

筛选蛋白质链的过程：
1.挑选出于其他链相似最少的，把这条链移到新的存放最终结果的链的列表中，从矩阵中剔除掉所有跟它相似度大于30%的蛋白质链以及这条链的信息；
2.新的矩阵继续1的操作，直到矩阵没有字符1为止，剩下全为0的矩阵中链表信息添加到最终结果的链的列表中
3.最后保存的链的列表信息，即为最终的非冗余数据集。
具体Python代码如下：

#对矩阵记性筛选def select(matrix):    #选出和最小的一行，记录ID    minid=0    minsum=sum(matrix[0,1:226])    for i in range(matrix.shape[0]):        su=sum(matrix[i,1:226])        if su<minsum:            minid=i            minsum=su    print minid,minsum    #选出保留的一个pdbchain之后，如果与这一个pdbchain相关的其他的蛋白质链都删除，即这个ID行，列中是1的对应的那一行和一列    dellist=[minid+1]    remain=matrix[minid,0]    for j in range(1,matrix.shape[1]):        if matrix[minid,j]==1:            dellist.append(j)    if len(dellist)>0:        dellistrow=[]        for k in dellist:            dellistrow.append(k-1)        try:            matrix=np.delete(matrix,dellist,axis=1)#删除列            matrix=np.delete(matrix, dellistrow, axis=0)        except:            print "delete error"    return matrix, remain

具体完整的代码，待我处理好所有的数据以及工作，会一同上传。这只是小白的一个实验步骤，有不对的地方希望大佬们指教，如对你有帮助，请点赞哟~~