胰蛋白酶消化（切割）无法使用正则表达式

Question

我正在尝试用Python编写一个程序，来模拟蛋白质序列的胰蛋白酶切割。胰蛋白酶的切割规则是：在R或K后面切割，但在P前面不切割。也就是说，胰蛋白酶会在每个K或R后面切割蛋白质序列，前提是这个K或R后面没有P。

举个例子，给定序列 MVPPPPSRGGAAKPGQLGRSLGPLLLLLRPEEPEDGDREICSESK，切割后应该得到这4个序列（肽段）：

MVPPPPSR
GGAAKPGQLGR
SLGPLLLLLRPEEPEDGDR
EICSESK 

注意，在第二个肽段中，K后面没有切割（因为K后面跟着P），在第三个肽段中，R后面也没有切割（同样因为R后面跟着P）。

我已经用Python写了这段代码，但效果不好。有没有什么方法可以更有效地实现这个正则表达式？

    # Open the file and read it line by line.

    myprotein = open(raw_input('Enter input filename: '),'r')
    if  os.path.exists("trypsin_digest.txt"):
        os.remove("trypsin_digest.txt")
    outfile = open("trypsin_digest.txt",'w+')

    for line in myprotein:
        protein = line.rstrip()
        protein = re.sub('(?<=[RK])(?=[^P])','', protein)

    for peptide in protein:
        outfile.write(peptide)
    print 'results written to:\n', os.getcwd() +'\ trypsin_digest.txt'

这是我让它正常工作的方式

   myprotein = open(raw_input('Enter input filename: '),'r')
   my_protein = []

   for protein in myprotein:
   myprotein = protein.rstrip('\n')
   my_protein.append(myprotein)
   my_pro = (''.join(my_protein))

   #cleaves sequence    
   peptides = re.sub(r'(?<=[RK])(?=[^P])','\n', my_pro)
   print peptides

蛋白质序列：

MVPPPPSRGGAAKPGQLGRSLGPLLLLLRPEEPEDGDREICSESK MVPPPPSRGGAAKPGQLGRSLGPLLLLLRPEEPEDGDREICSESK MVPPPPSRGGAAKPGQLGRSLGPLLLLLRPEEPEDGDREICSESK

输出（胰蛋白酶切割位置）或肽段

MVPPPPSR
GGAAKPGQLGR
SLGPLLLLLRPEEPEDGDR
EICSESK
MVPPPPSR
GGAAKPGQLGR
SLGPLLLLLRPEEPEDGDR
EICSESK
MVPPPPSR
GGAAKPGQLGR
SLGPLLLLLRPEEPEDGDR
EICSESK

Answer 1

我认为下面这个正则表达式可以满足你的需求：

([KR]?[^P].*?[KR](?!P))

下面是来自 pythonregexp 的结果

>>> regex = re.compile("([KR]?[^P].*?[KR](?!P))")
>>> r = regex.search(string)
>>> r
<_sre.SRE_Match object at 0xb1a9f49eb4111980>
>>> regex.match(string)
<_sre.SRE_Match object at 0xb1a9f49eb4102980>

# List the groups found
>>> r.groups()
(u'MVPPPPSR',)

# List the named dictionary objects found
>>> r.groupdict()
{}

# Run findall
>>> regex.findall(string)
[u'MVPPPPSR', u'GGAAKPGQLGR', u'SLGPLLLLLRPEEPEDGDR', u'EICSESK']

Answer 2

编辑 经过一些小调整，你的正则表达式效果很好：

在你的评论中，你提到你在一个文件里有多个序列（我们称这个文件为 sequences.dat）：

$ cat sequences.dat
MVPPPPSRGGAAKPGQLGRSLGPLLLLLRPEEPEDGDREICSESK
MVPPPPSRGGAAKPGQLGRSLGPLLLLLRPEEPEDGDREICSESK
MVPPPPSRGGAAKPGQLGRSLGPLLLLLRPEEPEDGDREICSESK

>>> with open('sequences.dat') as f:
    s = f.read()

>>> print(s)
MVPPPPSRGGAAKPGQLGRSLGPLLLLLRPEEPEDGDREICSESK
MVPPPPSRGGAAKPGQLGRSLGPLLLLLRPEEPEDGDREICSESK
MVPPPPSRGGAAKPGQLGRSLGPLLLLLRPEEPEDGDREICSESK

>>> protein = re.sub(r'(?<=[RK])(?=[^P])','\n', s, re.DOTALL)

>>> protein.split()
['MVPPPPSR', 'GGAAKPGQLGR', 'SLGPLLLLLRPEEPEDGDR', 'EICSESK', 'MVPPPPSR', 'GGAAKPGQLGR', 'SLGPLLLLLRPEEPEDGDR', 'EICSESK', 'MVPPPPSR', 'GGAAKPGQLGR', 'SLGPLLLLLRPEEPEDGDR', 'EICSESK']

>>> print protein
MVPPPPSR
GGAAKPGQLGR
SLGPLLLLLRPEEPEDGDR
EICSESK

MVPPPPSR
GGAAKPGQLGR
SLGPLLLLLRPEEPEDGDR
EICSESK

MVPPPPSR
GGAAKPGQLGR
SLGPLLLLLRPEEPEDGDR
EICSESK

Answer 3

正则表达式很不错，但这里有一个使用普通Python的解决方案。因为你是在寻找基础中的子序列，所以把这个做成一个生成器是很合理的，这样可以逐个输出片段。

example = 'MVPPPPSRGGAAKPGQLGRSLGPLLLLLRPEEPEDGDREICSESK'

def trypsin(bases):
    sub = ''
    while bases:
        k, r = bases.find('K'), bases.find('R')
        cut = min(k, r)+1 if k > 0 and r > 0 else max(k, r)+1
        sub += bases[:cut]
        bases = bases[cut:]
        if not bases or bases[0] != 'P':
            yield sub
            sub = ''


print list(trypsin(example))