DNA功能域是一段保守的有特定功能的一小段序列，在生物进化过程中保持保守，长度一般不超过520bp，通常和转录因子结合，以方便快速定位功能蛋白。DNA功能域能出现在基因的不同位置，出现的次数也不确定。
本软件使用的是分支界定算法，即将DNA序列变为树状结构，每个节点下面有4个分支，分别为A，T，C，G。需要使用者输入想发现的DNA功能域长度，然后计算hammingDistance，使用通常的贪婪算法需要遍历树的每个节点，本软件使用suffix即后缀算法，如"mississip"，其实就是在字典树中插入mississip,ississip，sissip,issip,ssip,sip,ip,p。DNA序列ATCGGTCG由suffix G，CG，TCG，GTCG，GGTCG，TCGGTCG，ATCGGTCG组成。计算hammingDistance评分的时候如果序列的后缀部分已经不好，那么就没有必要计算完整的部分，在树的结构里就会跳过即bypass。假如需要比对的有5条DNA序列，那么会将其用kmer构建树结构，代码如下。
import re,sys
def nextvertix(array,i,L,k): # 去树结构中的下一个节点
if i<L:
array.append(1)
return array,i+1
else:
for j in range(L,0,-1):
if array[j-1]<k:
array[j-1]+=1
return array,len(array)
else:
array[j-1]=1
array[-1:]=[]
return array,0
q = [2,4,4,4]
z = nextvertix(q,len(q),5,4)
#print(z)
def byPass(a,i,L,k): #直接跳到下个同等级的下一个节点
for j in range(i,0,-1):
if a[j-1]<k:
a[j-1]+=1
return a,len(a)
a[j-1]=1
a[-1:]=[]
return a,0
a = [2,4,4,3,3]
d = byPass(a,5,5,4)
#print(d)
def kmer(length,string):# 将DNA序列分割成kmer
lis = []
for i in range(0,len(string)-length+1):
lis.append(string[i:length+i])
return lis
def lettreDiffrent(l,n):#计算2条序列之间的不同字母数量
dis = 0
for i,j in zip(l,n):
if i!=j:
dis+=1
return dis
def hammingDistance(s,DNA):#计算hamming distance
liss = []
for i in DNA:
lis = kmer(len(s),i)
if s in lis:
liss.append(0)
else:
liss.append(min([lettreDiffrent(s, i) for i in lis]))
return sum(liss)
b = 'GCTTAGCCGGCCTGGCCAAT'
d = 'ACTTGCAGAGCGGCGGTTAA'
f = 'GGTTGCGTCCCAACTGGACC'
'''b = 'GCTTAGCCGGCCTGGCCAATGG'
d = 'AC'
f = 'GGTTGCGTCCCAACTGGACCAA'
i = 'AGTTGCGTCACTTCCTGGGCC'
'''
DNA = [b,d,f]
#DNA=['CGGGCCTG','ACCTGGCA','CACCTGGC','GCCAACGT']
def arrayToStr(array):#将数字变成字母
dic = {1:'A',2:'C',3:'G',4:'T'}
a = ''.join([dic[i] for i in array])
return a
def branchAndBound(DNA,l):# 主函数
liste = []
s = [1]
bestDis = 10000000
i = 1
bestWord= ''
while i>0:
if i<l:
prefix = s[:i]
#print(s,i)
suffix = bestWord[i:]
#print(arrayToStr(prefix),suffix,bestWord)
#print(prefix)
#print(arrayToStr(prefix))
optimalDis = hammingDistance(arrayToStr(prefix),DNA)
optimalSuffixDis = hammingDistance(suffix,DNA)
#print(optimalDis,bestDis)
if optimalDis+optimalSuffixDis>bestDis:
s,i= byPass(s,len(s),l,4)
#print('bypass')
#print(prefix)
#print(s)
#print('o')
else:
#print(prefix)
s,i= nextvertix(s,i,l,4)
#print(s)
else:
word = arrayToStr(s)
#print(s)
#print(word)
liste.append(word)
if hammingDistance(word,DNA)<bestDis:
bestDis=hammingDistance(word,DNA)
#print(bestDis)
bestWord = word
#print(bestWord)
#print(word)
s,i= nextvertix(s,i,l,4)
#print(i)
return bestWord, liste
import time
print('你好，欢迎使用由高端设计的DNA功能域搜索器，该程序使用的是分支界定的suffix算法')
print('请输入你想发现功能域的DNA序列，请每次输入一条DNA序列，按回车键继续下一条DNA序列输入')
print('='*37,'菜单','='*37)
print('1.输入DNA序列,输入q结束输入,接下来输入功能域长度')
print('2.关闭程序')
lis=[]
userinput=input('>>>>>')
if userinput=='2':
sys.exit(0)
while userinput!='2':
try:
userinput=input('请输入DNA序列或者输入q(Q)退出序列输入:').upper()
if userinput=='':
raise Exception
if userinput=='Q' and len(lis)>=2:
break
elif userinput=='Q' and len(lis)<2:
print('你好，请至少输入2条DNA序列，程序才能进行比对')
result=re.findall(r'[^ATCG]',userinput)
if result!=[] and result!='[Q]':
raise Exception
if userinput!='Q':
lis.append(userinput)
except Exception:
print('请输入有效的DNA序列')
print(lis)
min_len=min([len(i) for i in lis])
print(min_len)
loop=True
while loop:
try:
lengths=int(input('>>>>>长度:'))
if lengths>min_len:
raise Exception
loop=False
except Exception:
print('对不起，你希望的功能域长度大于你序列最短的长度，请重新输入')
print('开始运算<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<')
start_time = time.time()
a,liste = branchAndBound(lis,lengths)
print('发现功能域',a)
#print(len(liste))
print('寻找功能域消耗的运算时间为:',round(time.time()-start_time,3),'秒')
软件也可以使用python3 pip install suffixdg下载使用。
pypi.org个人项目链接：项目链接
