|
我正在研究一个简单的生物信息学问题.我有一个有效的解决方案,但这是非常低效的.我怎样才能提高效率？ 问题： 在字符串g中找到长度为k的模式,假设k-mer可能具有最多d个不匹配. 这些字符串和模式都是基因组的 – 所以我们的可能字符集是{A,T,C,G}. 我将调用函数FrequentWordsMismatch(g,k,d). 所以,这里有一些有用的例子： FrequentWordsMismatch(‘AAAAAAAAAA’,2,1)→[‘AA’,’CA’,’GA’,’TA’,’AC’,’AG’,’AT’] 这是一个更长的例子,如果你实现这个并想测试： FrequentWordsMisMatch( ‘CACAGTAGGCGCCGGCACACACAGCCCCGGGCCCCGGGCCGCCCCGGGCCGGCGGCCGCCGGCGCCGGCACACCGGCACAGCCGTACCGGCACAGTAGTACCGGCCGGCCGGCACACCGGCACACCGGGTACACACCGGGGCGCACACACAGGCGGGCGCCGGGCCCCGGGCCGTACCGGGCCGCCGGCGGCCCACAGGCGCCGGCACAGTACCGGCACACACAGTAGCCCACACACAGGCGGGCGGTAGCCGGCGCACACACACACAGTAGGCGCACAGCCGCCCACACACACCGGCCGGCCGGCACAGGCGGGCGGGCGCACACACACCGGCACAGTAGTAGGCGGCCGGCGCACAGCC’,10,2)→[ ‘GCACACAGAC’, ‘GCGCACACAC’] 凭借我天真的解决方案,第二个例子很容易花费大约60秒,尽管第一个例子非常快. 天真的解决方案： 我的想法是,对于g中的每个k长度段,找到每个可能的“邻居”(例如,具有多达d个不匹配的其他k长度段)并将这些邻居作为键添加到字典中.然后我计算每个邻居kmers出现在字符串g中的次数,并将其记录在字典中. 显然,这是一种有点糟糕的方式,因为邻居的数量像k和d一样疯狂地扩展,并且必须通过每个邻居扫描字符串使得这种实现非常缓慢.但是,这就是我要求帮助的原因. 我会把我的代码放在下面.打开包装肯定会有很多新手错误,所以感谢您的时间和精力. def FrequentWordsMismatch(g, k, d):
'''
Finds the most frequent k-mer patterns in the string g, given that those
patterns can mismatch amongst themselves up to d times
g (String): Collection of {A, T, C, G} characters
k (int): Length of desired pattern
d (int): Number of allowed mismatches
'''
counts = {}
answer = []
for i in range(len(g) - k 1):
kmer = g[i:i k]
for neighborkmer in Neighbors(kmer, d):
counts[neighborkmer] = Count(neighborkmer, g, d)
maxVal = max(counts.values())
for key in counts.keys():
if counts[key] == maxVal:
answer.append(key)
return(answer)
def Neighbors(pattern, d):
'''
Find all strings with at most d mismatches to the given pattern
pattern (String): Original pattern of characters
d (int): Number of allowed mismatches
'''
if d == 0:
return [pattern]
if len(pattern) == 1:
return ['A', 'C', 'G', 'T']
answer = []
suffixNeighbors = Neighbors(pattern[1:], d)
for text in suffixNeighbors:
if HammingDistance(pattern[1:], text) < d:
for n in ['A', 'C', 'G', 'T']:
answer.append(n text)
else:
answer.append(pattern[0] text)
return(answer)
def HammingDistance(p, q):
'''
Find the hamming distance between two strings
p (String): String to be compared to q
q (String): String to be compared to p
'''
ham = 0 abs(len(p)-len(q))
for i in range(min(len(p), len(q))):
if p[i] != q[i]:
ham = 1
return(ham)
def Count(pattern, g, d):
'''
Count the number of times that the pattern occurs in the string g,
allowing for up to d mismatches
pattern (String): Pattern of characters
g (String): String in which we're looking for pattern
d (int): Number of allowed mismatches
'''
return len(MatchWithMismatch(pattern, g, d))
def MatchWithMismatch(pattern, g, d):
'''
Find the indicies at which the pattern occurs in the string g,
allowing for up to d mismatches
pattern (String): Pattern of characters
g (String): String in which we're looking for pattern
d (int): Number of allowed mismatches
'''
answer = []
for i in range(len(g) - len(pattern) 1):
if(HammingDistance(g[i:i len(pattern)], pattern) <= d):
answer.append(i)
return(answer)
更多测试 FrequentWordsMismatch('ACGTTGCATGTCGCATGATGCATGAGAGCT', 4, 1) → ['ATGC', 'ATGT', 'GATG']
FrequentWordsMismatch('AGTCAGTC', 4, 2) → ['TCTC', 'CGGC', 'AAGC', 'TGTG', 'GGCC', 'AGGT', 'ATCC', 'ACTG', 'ACAC', 'AGAG', 'ATTA', 'TGAC', 'AATT', 'CGTT', 'GTTC', 'GGTA', 'AGCA', 'CATC']
FrequentWordsMismatch('AATTAATTGGTAGGTAGGTA', 4, 0) → ["GGTA"]
FrequentWordsMismatch('ATA', 3, 1) → ['GTA', 'ACA', 'AAA', 'ATC', 'ATA', 'AGA', 'ATT', 'CTA', 'TTA', 'ATG']
FrequentWordsMismatch('AAT', 3, 0) → ['AAT']
FrequentWordsMismatch('TAGCG', 2, 1) → ['GG', 'TG']
解决方法: 问题描述在几个方面是模棱两可的,所以我将通过这些例子.您似乎想要字母表中的所有k长度字符串(A,C,G,T),使得与g的连续子串的匹配数最大 – 其中“匹配”表示最多与字符相同d字符不等式. 我忽略了你的HammingDistance()函数即使在输入长度不同的情况下也会产生一些东西,主要是因为它对我没有多大意义;-),但部分是因为不需要得到你想要的结果你提供的任何例子. 下面的代码生成了所有示例中所需的结果,从而产生了您给出的输出列表的排列.如果你想要规范输出,我建议在返回之前对输出列表进行排序. 该算法非常简单,但依赖于itertools来“以C速度”进行重组合提升.所有的例子都在第二个总数下运行良好. 对于g的每个长度-k连续子串,考虑d个不同索引位置的所有组合(k,d)组.有四种方法用来自{A,C,G,T}的字母填充那些索引位置,并且每种方式都是“一种模式”,其匹配具有最多d个差异的子字符串.通过记住已生成的模式来清除重复项;这比制作英雄的努力更快地生成只有独特的模式开始. 所以,总的来说,时间要求是O(len(g)* k ** d * 4 ** d)= O(len(g)*(4 * k)** d,其中k ** d是,对于k和d的相当小的值,二项式系数组合(k,d)被夸大了.重要的是要注意的是 – 毫不奇怪 – 它在d中是指数的. def fwm(g, k, d):
from itertools import product, combinations
from collections import defaultdict
all_subs = list(product("ACGT", repeat=d))
all_ixs = list(combinations(range(k), d))
patcount = defaultdict(int)
for starti in range(len(g)):
base = g[starti : starti k]
if len(base) < k:
break
patcount[base] = 1
seen = set([base])
basea = list(base)
for ixs in all_ixs:
saved = [basea[i] for i in ixs]
for newchars in all_subs:
for i, newchar in zip(ixs, newchars):
basea[i] = newchar
candidate = "".join(basea)
if candidate not in seen:
seen.add(candidate)
patcount[candidate] = 1
for i, ch in zip(ixs, saved):
basea[i] = ch
maxcount = max(patcount.values())
return [p for p, c in patcount.items() if c == maxcount]
编辑：独特地生成模式 而不是通过保留迄今为止看到的那些副本来清除重复项,它足以直接防止生成重复项.事实上,下面的代码更短更简单,虽然有点微妙.作为少量冗余工作的回报,有一些对inner()函数的递归调用.哪种方式更快似乎取决于具体的输入. def fwm(g, k, d):
from collections import defaultdict
patcount = defaultdict(int)
alphabet = "ACGT"
allbut = {ch: tuple(c for c in alphabet if c != ch)
for ch in alphabet}
def inner(i, rd):
if not rd or i == k:
patcount["".join(base)] = 1
return
inner(i 1, rd)
orig = base[i]
for base[i] in allbut[orig]:
inner(i 1, rd-1)
base[i] = orig
for i in range(len(g) - k 1):
base = list(g[i : i k])
inner(0, d)
maxcount = max(patcount.values())
return [p for p, c in patcount.items() if c == maxcount]
来源：http://www./content-1-202701.html
|
