DNA 字符串中随机基序搜索的意外输出
Unexpected output in randomized motif search in DNA strings
我有以下 t=5 DNA 字符串：
DNA = '''CGCCCCTCTCGGGGGTGTTCAGTAAACGGCCA
GGGCGAGGTATGTGTAAGTGCCAAGGTGCCAG
TAGTACCGAGACCGAAAGAAGTATACAGGCGT
TAGATCAAGTTTCAGGTGCACGTCGGTGAACC
AATCCACCAGCTCCACGTGCAATGTTGGCCTA'''
k = 8
t = 5
我正在尝试使用拉普拉斯变换从字符串集合中找到长度为 k=8 的最佳图案，以从每个 t 字符串中随机抽取长度为 k 的块。
我的辅助函数如下：
def window(s, k):
for i in range(1 + len(s) - k):
yield s[i:i+k]
def HammingDistance(seq1, seq2):
if len(seq1) != len(seq2):
raise ValueError('Undefined for sequences of unequal length.')
return sum(ch1 != ch2 for ch1, ch2 in zip(seq1, seq2))
def score(motifs):
score = 0
for i in range(len(motifs[0])):
motif = ''.join([motifs[j][i] for j in range(len(motifs))])
score += min([HammingDistance(motif, homogeneous*len(motif)) for homogeneous in 'ACGT'])
return score
def profile(motifs):
prof = []
for i in range(len(motifs[0])):
col = ''.join([motifs[j][i] for j in range(len(motifs))])
prof.append([float(col.count(nuc))/float(len(col)) for nuc in 'ACGT'])
return prof
def profile_most_probable_kmer(dna, k, prof):
dna = dna.splitlines()
nuc_loc = {nucleotide:index for index,nucleotide in enumerate('ACGT')}
motif_matrix = []
max_prob = [-1, None]
for i in range(len(dna)):
motif_matrix.append(max_prob)
for i in range(len(dna)):
for chunk in window(dna[i],K):
current_prob = 1
for j, nuc in enumerate(chunk):
current_prob*=prof[j][nuc_loc[nuc]]
if current_prob>motif_matrix[i][0]:
motif_matrix[i] = [current_prob,chunk]
return list(list(zip(*motif_matrix))[1])
def profile_with_pseudocounts(motifs):
prof = []
for i in range(len(motifs[0])):
col = ''.join([motifs[j][i] for j in range(len(motifs))])
prof.append([float(col.count(nuc)+1)/float(len(col)+4) for nuc in 'ACGT'])
return prof
from random import randint
def SampleMotifs(Dna,k,t):
Dna = Dna.splitlines()
BestMotifs = []
for line in Dna:
position = randint(0,len(line)-k)
BestMotifs.append(line[position:position+k])
return BestMotifs
def motifs_from_profile(profile, dna, k):
return [profile_most_probable_kmer(seq,k,profile) for seq in dna]
def randomized_motif_search(dna,k,t):
from random import randint
dna = dna.splitlines()
rand_ints = [randint(0,len(dna[0])-k)) for a in range(len(dna))]
motifs = [dna[i][r:r+k] for i,r in enumerate(rand_ints)]
best_score = [score(motifs), motifs]
while True:
current_profile = profile_with_pseudocounts(motifs)
motifs = motifs_from_profile(current_profile, dna, k)
current_score = score(motifs)
if current_score < best_score[0]:
best_score = [current_score, motifs]
else:
return best_score[1]
def Laplace(dna,k,t):
i = 0
LastMotifs = randomized_motif_search(dna,k,t)
while i < 1000:
try:
BestMotifs = randomized_motif_search(dna,k,t)
if score(BestMotifs)<score(LastMotifs):
LastMotifs = BestMotifs
except:
pass
i+=1
print(*LastMotifs)
我应该得到的输出是：
TCTCGGGG
CCAAGGTG
TACAGGCG
TTCAGGTG
TCCACGTG
我每次都得到不同的输出，这与我期望的带有随机元素的方法一样，但是当我迭代 1000 次时它应该会收敛，并且只在分数较低时更新我最好的图案。我不得不在 laplace 中放入一个错误处理程序，因为我在调用 randomized_motif_search(dna,k,t) 时获得了一个索引，这告诉我这可能是问题的根源。在过去的两天里，我一直在搜索代码以确保一切都是正确的，但事实上我得到了错误的答案或以下错误：
---------------------------------------------------------------------------
IndexError Traceback (most recent call last)
<ipython-input-811-ee735449bb8e> in <module>
----> 1 Laplace(DNA,K,T)
<ipython-input-810-31301d79eb95> in Laplace(dna, k, t)
3 LastMotifs = randomized_motif_search(dna,k,t)
4 while i < 2000:
----> 5 BestMotifs = randomized_motif_search(dna,k,t)
6 if score(BestMotifs)<score(LastMotifs):
7 LastMotifs = BestMotifs
<ipython-input-809-43600d882734> in randomized_motif_search(dna, k, t)
8 while True:
9 current_profile = profile_with_pseudocounts(motifs)
---> 10 motifs = motifs_from_profile(current_profile, dna, k)
11 current_score = score(motifs)
12 if current_score < best_score[0]:
<ipython-input-408-7c866045d839> in motifs_from_profile(profile, dna, k)
1 def motifs_from_profile(profile, dna, k):
----> 2 return [profile_most_probable_kmer(seq,k,profile) for seq in dna]
<ipython-input-408-7c866045d839> in <listcomp>(.0)
1 def motifs_from_profile(profile, dna, k):
----> 2 return [profile_most_probable_kmer(seq,k,profile) for seq in dna]
<ipython-input-795-56f83ba5ee2b> in profile_most_probable_kmer(dna, k, prof)
25 current_prob = 1
26 for j, nuc in enumerate(chunk):
---> 27 current_prob*=prof[j][nuc_loc[nuc]]
28 if current_prob>motif_matrix[i][0]:
29 motif_matrix[i] = [current_prob,chunk]
IndexError: list index out of range
这有点令人烦恼。实际帮助将不胜感激。
编辑：问题是我对从 dna 字符串中抽取随机图案的随机整数进行索引，而我的 motifs_from_profile 函数 return 是一个列表列表，而不仅仅是一个列表代码有点依赖。我已经更新了以下功能：虽然这些修复解决了代码中引发错误的问题，但我现在每次 运行 Laplace 函数时都会得到一个输出，但结果不是即使我在第一次迭代时输入正确答案，我的期望也是如此。我将尽我最大的努力调试评分函数中发生的事情，并回顾我猜想的文献。也许来自社区的一些更模糊的意见会有所帮助，但谁知道呢？
更新后的 randomized_motif_search 是：
def randomized_motif_search(dna,k,t):
from random import randint
from itertools import chain
dna = dna.splitlines()
# Randomly generate k-mers from each sequence in the dna list.
rand_ints = [randint(0,len(dna[0])-(k)) for a in range(len(dna))]
motifs = [dna[i][r:r+k] for i,r in enumerate(rand_ints)]
best_score = [score(motifs), motifs]
while True:
current_profile = profile_with_pseudocounts(motifs)
mfp = motifs_from_profile(current_profile,dna,k)
motifs = []
for i in range(len(mfp)):
motifs.append(mfp[i][0])
current_score = score(motifs)
if current_score < best_score[0]:
best_score = [current_score, motifs]
else:
return best_score[1]
和新的 Laplace:
def Laplace(dna,k,t):
i = 0
LastMotifs = randomized_motif_search(dna,k,t)
while i < 1000:
BestMotifs = randomized_motif_search(dna,k,t)
if score(BestMotifs)<score(LastMotifs):
LastMotifs = BestMotifs
i+=1
print(*LastMotifs)
编辑：亲爱的杂志，
我弄乱了得分方法并想出了如何 return 主题得分，但仍然没有得到问题的正确答案。我正式卡在这里是更新后的 score 具有正确索引的函数代码：
def score(motifs):
score = 0
for i in range(len(motifs[0])):
motif = ''.join([Motifs[j][i] for j in range(len(Motifs))])
score+=min([HammingDistance(motif,homogenous*len(motif)) for homogenous in 'ACGT'])
return score
我想通了！所有的索引都源于对 splitlines() 的调用，我在许多辅助函数的开头进行了调用。我也没有完全调试我的 score() 函数，所以我的评分并没有按照文献要求的方式进行。
所有辅助函数如下：
def window(s, k):
for i in range(1 + len(s) - k):
yield s[i:i+k]
def HammingDistance(seq1, seq2):
if len(seq1) != len(seq2):
raise ValueError('Undefined for sequences of unequal length.')
return sum(ch1 != ch2 for ch1, ch2 in zip(seq1, seq2))
def score(motifs):
score = 0
for i in range(len(motifs[0])):
motif = ''.join(motifs[j][i] for j in range(len(motifs)))
score += min([HammingDistance(motif,homogenous*len(motif)) for homogenous in 'ACGT'])
return score
def profile(motifs):
prof = []
for i in range(len(motifs[0])):
col = ''.join([motifs[j][i] for j in range(len(motifs))])
prof.append([float(col.count(nuc))/float(len(col)) for nuc in 'ACGT'])
return prof
def profile_most_probable_kmer(dna, k, prof):
nuc_loc = {nucleotide:index for index,nucleotide in enumerate('ACGT')}
motif_matrix = []
max_prob = [-1, None]
for i in range(len(dna)):
motif_matrix.append(max_prob)
for i in range(len(dna)):
for chunk in window(dna[i],k):
current_prob = 1
for j, nuc in enumerate(chunk):
current_prob*=prof[j][nuc_loc[nuc]]
if current_prob>motif_matrix[i][0]:
motif_matrix[i] = [current_prob,chunk]
return list(list(zip(*motif_matrix))[1])
def profile_with_pseudocounts(motifs):
prof = []
for i in range(len(motifs[0])):
col = ''.join([motifs[j][i] for j in range(len(motifs))])
prof.append([float(col.count(nuc)+1)/float(len(col)+4) for nuc in 'ACGT'])
return prof
from random import randint
def SampleMotifs(Dna,k,t):
Dna = Dna.splitlines()
BestMotifs = []
for line in Dna:
position = randint(0,len(line)-k)
BestMotifs.append(line[position:position+k])
return BestMotifs
def randomized_motif_search(dna,k,t):
from random import randint
from itertools import chain
dna = dna.splitlines()
# Randomly generate k-mers from each sequence in the dna list.
rand_ints = [randint(0,len(dna[0])-(k)) for a in range(len(dna))]
motifs = [dna[i][r:r+k] for i,r in enumerate(rand_ints)]
best_score = [score(motifs), motifs]
while True:
current_profile = profile_with_pseudocounts(motifs)
motifs = profile_most_probable_kmer(dna,k,current_profile)
# motifs = []
# for i in range(len(mfp)):
# motifs.append(mfp[i][0])
current_score = score(motifs)
if current_score < best_score[0]:
best_score = [current_score, motifs]
else:
return best_score[1]
def Laplace(dna,k,t):
import random
random.seed(0)
i = 0
LastMotifs = randomized_motif_search(dna,k,t)
while i < 1000:
BestMotifs = randomized_motif_search(dna,k,t)
if score(BestMotifs)<score(LastMotifs):
LastMotifs = BestMotifs
i+=1
print('\n'.join(LastMotifs))
我有以下 t=5 DNA 字符串：
DNA = '''CGCCCCTCTCGGGGGTGTTCAGTAAACGGCCA
GGGCGAGGTATGTGTAAGTGCCAAGGTGCCAG
TAGTACCGAGACCGAAAGAAGTATACAGGCGT
TAGATCAAGTTTCAGGTGCACGTCGGTGAACC
AATCCACCAGCTCCACGTGCAATGTTGGCCTA'''
k = 8
t = 5
我正在尝试使用拉普拉斯变换从字符串集合中找到长度为 k=8 的最佳图案，以从每个 t 字符串中随机抽取长度为 k 的块。
我的辅助函数如下：
def window(s, k):
for i in range(1 + len(s) - k):
yield s[i:i+k]
def HammingDistance(seq1, seq2):
if len(seq1) != len(seq2):
raise ValueError('Undefined for sequences of unequal length.')
return sum(ch1 != ch2 for ch1, ch2 in zip(seq1, seq2))
def score(motifs):
score = 0
for i in range(len(motifs[0])):
motif = ''.join([motifs[j][i] for j in range(len(motifs))])
score += min([HammingDistance(motif, homogeneous*len(motif)) for homogeneous in 'ACGT'])
return score
def profile(motifs):
prof = []
for i in range(len(motifs[0])):
col = ''.join([motifs[j][i] for j in range(len(motifs))])
prof.append([float(col.count(nuc))/float(len(col)) for nuc in 'ACGT'])
return prof
def profile_most_probable_kmer(dna, k, prof):
dna = dna.splitlines()
nuc_loc = {nucleotide:index for index,nucleotide in enumerate('ACGT')}
motif_matrix = []
max_prob = [-1, None]
for i in range(len(dna)):
motif_matrix.append(max_prob)
for i in range(len(dna)):
for chunk in window(dna[i],K):
current_prob = 1
for j, nuc in enumerate(chunk):
current_prob*=prof[j][nuc_loc[nuc]]
if current_prob>motif_matrix[i][0]:
motif_matrix[i] = [current_prob,chunk]
return list(list(zip(*motif_matrix))[1])
def profile_with_pseudocounts(motifs):
prof = []
for i in range(len(motifs[0])):
col = ''.join([motifs[j][i] for j in range(len(motifs))])
prof.append([float(col.count(nuc)+1)/float(len(col)+4) for nuc in 'ACGT'])
return prof
from random import randint
def SampleMotifs(Dna,k,t):
Dna = Dna.splitlines()
BestMotifs = []
for line in Dna:
position = randint(0,len(line)-k)
BestMotifs.append(line[position:position+k])
return BestMotifs
def motifs_from_profile(profile, dna, k):
return [profile_most_probable_kmer(seq,k,profile) for seq in dna]
def randomized_motif_search(dna,k,t):
from random import randint
dna = dna.splitlines()
rand_ints = [randint(0,len(dna[0])-k)) for a in range(len(dna))]
motifs = [dna[i][r:r+k] for i,r in enumerate(rand_ints)]
best_score = [score(motifs), motifs]
while True:
current_profile = profile_with_pseudocounts(motifs)
motifs = motifs_from_profile(current_profile, dna, k)
current_score = score(motifs)
if current_score < best_score[0]:
best_score = [current_score, motifs]
else:
return best_score[1]
def Laplace(dna,k,t):
i = 0
LastMotifs = randomized_motif_search(dna,k,t)
while i < 1000:
try:
BestMotifs = randomized_motif_search(dna,k,t)
if score(BestMotifs)<score(LastMotifs):
LastMotifs = BestMotifs
except:
pass
i+=1
print(*LastMotifs)
我应该得到的输出是：
TCTCGGGG
CCAAGGTG
TACAGGCG
TTCAGGTG
TCCACGTG
我每次都得到不同的输出，这与我期望的带有随机元素的方法一样，但是当我迭代 1000 次时它应该会收敛，并且只在分数较低时更新我最好的图案。我不得不在 laplace 中放入一个错误处理程序，因为我在调用 randomized_motif_search(dna,k,t) 时获得了一个索引，这告诉我这可能是问题的根源。在过去的两天里，我一直在搜索代码以确保一切都是正确的，但事实上我得到了错误的答案或以下错误：
---------------------------------------------------------------------------
IndexError Traceback (most recent call last)
<ipython-input-811-ee735449bb8e> in <module>
----> 1 Laplace(DNA,K,T)
<ipython-input-810-31301d79eb95> in Laplace(dna, k, t)
3 LastMotifs = randomized_motif_search(dna,k,t)
4 while i < 2000:
----> 5 BestMotifs = randomized_motif_search(dna,k,t)
6 if score(BestMotifs)<score(LastMotifs):
7 LastMotifs = BestMotifs
<ipython-input-809-43600d882734> in randomized_motif_search(dna, k, t)
8 while True:
9 current_profile = profile_with_pseudocounts(motifs)
---> 10 motifs = motifs_from_profile(current_profile, dna, k)
11 current_score = score(motifs)
12 if current_score < best_score[0]:
<ipython-input-408-7c866045d839> in motifs_from_profile(profile, dna, k)
1 def motifs_from_profile(profile, dna, k):
----> 2 return [profile_most_probable_kmer(seq,k,profile) for seq in dna]
<ipython-input-408-7c866045d839> in <listcomp>(.0)
1 def motifs_from_profile(profile, dna, k):
----> 2 return [profile_most_probable_kmer(seq,k,profile) for seq in dna]
<ipython-input-795-56f83ba5ee2b> in profile_most_probable_kmer(dna, k, prof)
25 current_prob = 1
26 for j, nuc in enumerate(chunk):
---> 27 current_prob*=prof[j][nuc_loc[nuc]]
28 if current_prob>motif_matrix[i][0]:
29 motif_matrix[i] = [current_prob,chunk]
IndexError: list index out of range
这有点令人烦恼。实际帮助将不胜感激。
编辑：问题是我对从 dna 字符串中抽取随机图案的随机整数进行索引，而我的 motifs_from_profile 函数 return 是一个列表列表，而不仅仅是一个列表代码有点依赖。我已经更新了以下功能：虽然这些修复解决了代码中引发错误的问题，但我现在每次 运行 Laplace 函数时都会得到一个输出，但结果不是即使我在第一次迭代时输入正确答案，我的期望也是如此。我将尽我最大的努力调试评分函数中发生的事情，并回顾我猜想的文献。也许来自社区的一些更模糊的意见会有所帮助，但谁知道呢？
更新后的 randomized_motif_search 是：
def randomized_motif_search(dna,k,t):
from random import randint
from itertools import chain
dna = dna.splitlines()
# Randomly generate k-mers from each sequence in the dna list.
rand_ints = [randint(0,len(dna[0])-(k)) for a in range(len(dna))]
motifs = [dna[i][r:r+k] for i,r in enumerate(rand_ints)]
best_score = [score(motifs), motifs]
while True:
current_profile = profile_with_pseudocounts(motifs)
mfp = motifs_from_profile(current_profile,dna,k)
motifs = []
for i in range(len(mfp)):
motifs.append(mfp[i][0])
current_score = score(motifs)
if current_score < best_score[0]:
best_score = [current_score, motifs]
else:
return best_score[1]
和新的 Laplace:
def Laplace(dna,k,t):
i = 0
LastMotifs = randomized_motif_search(dna,k,t)
while i < 1000:
BestMotifs = randomized_motif_search(dna,k,t)
if score(BestMotifs)<score(LastMotifs):
LastMotifs = BestMotifs
i+=1
print(*LastMotifs)
编辑：亲爱的杂志，
我弄乱了得分方法并想出了如何 return 主题得分，但仍然没有得到问题的正确答案。我正式卡在这里是更新后的 score 具有正确索引的函数代码：
def score(motifs):
score = 0
for i in range(len(motifs[0])):
motif = ''.join([Motifs[j][i] for j in range(len(Motifs))])
score+=min([HammingDistance(motif,homogenous*len(motif)) for homogenous in 'ACGT'])
return score
我想通了！所有的索引都源于对 splitlines() 的调用，我在许多辅助函数的开头进行了调用。我也没有完全调试我的 score() 函数，所以我的评分并没有按照文献要求的方式进行。
所有辅助函数如下：
def window(s, k):
for i in range(1 + len(s) - k):
yield s[i:i+k]
def HammingDistance(seq1, seq2):
if len(seq1) != len(seq2):
raise ValueError('Undefined for sequences of unequal length.')
return sum(ch1 != ch2 for ch1, ch2 in zip(seq1, seq2))
def score(motifs):
score = 0
for i in range(len(motifs[0])):
motif = ''.join(motifs[j][i] for j in range(len(motifs)))
score += min([HammingDistance(motif,homogenous*len(motif)) for homogenous in 'ACGT'])
return score
def profile(motifs):
prof = []
for i in range(len(motifs[0])):
col = ''.join([motifs[j][i] for j in range(len(motifs))])
prof.append([float(col.count(nuc))/float(len(col)) for nuc in 'ACGT'])
return prof
def profile_most_probable_kmer(dna, k, prof):
nuc_loc = {nucleotide:index for index,nucleotide in enumerate('ACGT')}
motif_matrix = []
max_prob = [-1, None]
for i in range(len(dna)):
motif_matrix.append(max_prob)
for i in range(len(dna)):
for chunk in window(dna[i],k):
current_prob = 1
for j, nuc in enumerate(chunk):
current_prob*=prof[j][nuc_loc[nuc]]
if current_prob>motif_matrix[i][0]:
motif_matrix[i] = [current_prob,chunk]
return list(list(zip(*motif_matrix))[1])
def profile_with_pseudocounts(motifs):
prof = []
for i in range(len(motifs[0])):
col = ''.join([motifs[j][i] for j in range(len(motifs))])
prof.append([float(col.count(nuc)+1)/float(len(col)+4) for nuc in 'ACGT'])
return prof
from random import randint
def SampleMotifs(Dna,k,t):
Dna = Dna.splitlines()
BestMotifs = []
for line in Dna:
position = randint(0,len(line)-k)
BestMotifs.append(line[position:position+k])
return BestMotifs
def randomized_motif_search(dna,k,t):
from random import randint
from itertools import chain
dna = dna.splitlines()
# Randomly generate k-mers from each sequence in the dna list.
rand_ints = [randint(0,len(dna[0])-(k)) for a in range(len(dna))]
motifs = [dna[i][r:r+k] for i,r in enumerate(rand_ints)]
best_score = [score(motifs), motifs]
while True:
current_profile = profile_with_pseudocounts(motifs)
motifs = profile_most_probable_kmer(dna,k,current_profile)
# motifs = []
# for i in range(len(mfp)):
# motifs.append(mfp[i][0])
current_score = score(motifs)
if current_score < best_score[0]:
best_score = [current_score, motifs]
else:
return best_score[1]
def Laplace(dna,k,t):
import random
random.seed(0)
i = 0
LastMotifs = randomized_motif_search(dna,k,t)
while i < 1000:
BestMotifs = randomized_motif_search(dna,k,t)
if score(BestMotifs)<score(LastMotifs):
LastMotifs = BestMotifs
i+=1
print('\n'.join(LastMotifs))