如何用python实现蛋白质磷酸化位点的搜索

小白最近做磷酸化实验，进行蛋白多序列比对的时候发现NCBI BLAST中多序列联配不是很准确，于是参考python，编写了一个找特定磷酸化位点的程序。

#!coding:utf-8
import re    #导入正则表达式模块
motif = input('Please input motif:')  #输入磷酸化位点
regexp = re.compile(motif)  #运用正则表达式搜索motif
seq = input("Please input protein sequence:")  #输入蛋白序列
all = regexp.findall(seq)  #寻找所有的匹配
iter = regexp.finditer(seq)  #找到所有匹配并返回迭代器形式
list = []  #建立空列表
for s in iter:  #用for语句输出
    num = (s.start()+1) #由于s.start（）从0开始说以真正位点要加1
    print(str(s.start()+1)+' '+s.group())  #输出磷酸化位点位置和蛋白序列
    list.append(num)  #把所有的num 位置追加到list列表
#上面得到可以放到def模块里

if list == []:  #如果不匹配，list[]为空，输出No match
    print("No match!")
else:             #
    if len(list)== 1: #前两个比较特殊，不能用循环，否则会超出范围
        seq1 = seq[:list[0]-1]+'\033[1;31;40m'+ seq[list[0]-1:list[0]+1]+'\033[0m'+ seq[list[0]+1:]  #输入完整序列，磷酸化位点用红色字体打印
    if len(list) == 2:
        seq1 = seq[:list[0]-1]+'\033[1;31;40m'+ seq[list[0]-1:list[0]+1]+'\033[0m'+seq[list[0]+1:list[1]-1]+'\033[1;31;40m'+
seq[list[1]-1:list[1]+1]+'\033[0m'+ seq[list[1]+1:]


    if len(list) > 2:  #2以上可以用循环
        seq1 = seq[:list[0]-1]+'\033[1;31;40m'+ seq[list[0]-1:list[0]+1]+'\033[0m'+ seq[list[0]+1:list[1]-1]  #第一个序列为共用
        n=0    #设置一个计数器
        while n < len(list)-2: #如果n小于磷酸化位点个数-2

            n = n+1 #n不断加1
            seq1 = seq1+ '\033[1;31;40m'+seq[list[n]-1:list[n]+1]+ '\033[0m'+  seq[list[n]+1:list[n+1]-1]  #list从1开始，不断循环，位点两个氨基酸加位点和下一个位点的蛋白序列
            if n == (len(list)-2):  #直到n=磷酸化位点个数-2，停止循环。
                seq1 = seq1 +'\033[1;31;40m'+ seq[list[n+1]-1:list[n+1]+1]+'\033[0m'+seq[list[n+1]+1:] #加上最后磷酸化位点和最后一段序列
    print(seq1)  #输出所有标记序列

#至此，小程序已经完成，我们来运行一下试试