Python计算限制性核酸内切酶切割后的核酸片段及其片段分子量

       在基因工程中的重组DNA构建的实验中，往往需要用限制性核酸内切酶对DNA分子进行切割，已得到目的基因，但DNA分子上可能含有多个酶切位点，因此需要认为的计算酶切后得到的DNA片段（fragements）大小的分子量，以便后续电泳中可以根据Marker得到目的基因的片段，但对于酶切后得到的片段较多时，人工的计算片段的分子量较为耗时，因此这里特异根据《生物科学家的Python指南》一书，写了这篇文章以计算限制性核酸内切酶切割DNA分子后得到的片段分子量。

1.如何计算DNA片段的分子量？

这里用for循环写了一个函数，然后用一段DNA序列（

gcgatgctaggatccgc

）演示了函数的用法，具体代码如下：

#定义计算DNA片段的分子量大小的函数
def oligoMolecularWeight(seqence):
    dnaMolecularWeight = {'a': 313.2, 'c': 289.2, 't': 304.2, 'g': 329.2}
    molecularWeight = 0.0
    for base in seqence:
        molecularWeight += dnaMolecularWeight[base]
    return molecularWeight

#示例
print(oligoMolecularWeight('gcgatgctaggatccgc'))

 最后计算得到的结果为5273.399999999999。

2.find函数查找酶切位点在DNA上的位置，下面列了find函数运行后的结果：

mySequence = 'gcgatgctaggatccgcgatcgcgtacgatcgtacgcggtacggacggatccttctc'
print(mySequence.find('ggatcc'))
#运行结果：9
print(mySequence.find('ggatttt'))
#运行结果-1

由示例可以看出find函数在搜索字符串（‘ggatttt’） 时，因序列中无ggatttt，因此放回的是-1，表示未找到子字符串。

PS：find函数不能对“没有发现”返回零，因为这实际上是Python字符串中第一个位置的索引，如果返回零，则表示在该位置找到了这个子字符串。

3.while语句遍历整个DNA序列

mySequence = 'gcgatgctaggatccgcgatcgcgtacgatcgtacgcggtacggacggatccttctc'
found = 0
searchFrom = found
while found != -1:
    found = mySequence.find('ggatcc',searchFrom)
    if found != -1:
        print('Substring found at:',found)
    searchFrom = found +1


=====================================================================
运行结果：
Substring found at: 9
Substring found at: 46

4.合并计算分子量和酶切得到DNA片段的函数 

#创建限制性核酸内切酶的切割位点
restrictionEnzymes = {}
restrictionEnzymes['bamH1'] = ['ggatcc',0]
restrictionEnzymes['sma1'] = ['cccggg',2]

#定义计算DNA片段的分子量大小的函数
def oligoMolecularWeight(seqence):
    dnaMolecularWeight = {'a': 313.2, 'c': 289.2, 't': 304.2, 'g': 329.2}
    molecularWeight = 0.0
    for base in seqence:
        molecularWeight += dnaMolecularWeight[base]
    return molecularWeight

#定义酶切函数
def restrictionDigest(sequence,enzyme):
    #将限制性核酸内切酶的切割序列和切割位点分别赋值给motif和cutPosition
    motif = restrictionEnzymes[enzyme][0]
    cutPosition = restrictionEnzymes[enzyme][1]
    #提前创建用来存储切割后得到的酶切片段序列及其分子量的列表
    fragments = []
    #found表示上次找到的限制性核酸内切酶位点序列的起始位置
    found = 0
    #lastCut用来存储上次切割DNA的位置
    lastCut = found
    #searchFrom用来告诉我们下次从哪开始找出酶切序列
    searchFrom = lastCut
    #while循环遍历DNA序列
    while found != -1:
        found = sequence.find(motif,searchFrom)
        if found != -1:
            fragment = sequence[lastCut:found+cutPosition]
            wnt = oligoMolecularWeight(fragment)
            fragments.append((fragment,wnt))
        else:
            fragment = sequence[lastCut:]
            wnt = oligoMolecularWeight(fragment)
            fragments.append((fragment,wnt))
        lastCut = found+cutPosition
        searchFrom = lastCut+1
    return fragments


#用法示例
digestSquence = 'gcgatgctaggatccgcgatcgcgtacgatcgtacgcggtacggacggatccttctc'
print(restrictionDigest(digestSquence,'bamH1'))

     总结：大家只需要根据自己的需求添加所需的限制性核酸内切酶到字典，然后根据自己的序列进行运行就可以得到结果。