用Gibbs采样器进行主题搜索

Question

我是编程和生物信息学的初学者。所以，我会很感激你的理解。我尝试开发一个python脚本，使用Gibbs抽样进行主题搜索，如Coursera课程“在DNA中查找隐藏的消息”中所解释的。在使用过程中所提供的伪代码是：

GIBBSSAMPLER(Dna, k, t, N) 
    randomly select k-mers Motifs = (Motif1, …, Motift) in each string 
     from Dna 
    BestMotifs ← Motifs 
    for j ← 1 to N 
     i ← Random(t) 
     Profile ← profile matrix constructed from all strings in Motifs 
        except for Motifi 
     Motifi ← Profile-randomly generated k-mer in the i-th sequence 
     if Score(Motifs) < Score(BestMotifs) 
      BestMotifs ← Motifs 
    return BestMotifs 

问题描述：

CODE挑战：实现GIBBSSAMPLER。

输入：整数k，t和N，然后是字符串Dna的集合。 输出：从随机启动运行GIBBSSAMPLER（Dna，k，t，N）产生的字符串BestMotifs。请记住使用伪计数！

样品输入：

8 5 100 
CGCCCCTCTCGGGGGTGTTCAGTAACCGGCCA 
GGGCGAGGTATGTGTAAGTGCCAAGGTGCCAG 
TAGTACCGAGACCGAAAGAAGTATACAGGCGT 
TAGATCAAGTTTCAGGTGCACGTCGGTGAACC 
AATCCACCAGCTCCACGTGCAATGTTGGCCTA 

样本输出：

TCTCGGGG 
CCAAGGTG 
TACAGGCG 
TTCAGGTG 
TCCACGTG 

我跟着伪代码，以我的知识。这里是我的代码：

def BuildProfileMatrix(dnamatrix): 
    ProfileMatrix = [[1 for x in xrange(len(dnamatrix[0]))] for x in xrange(4)] 
    indices = {'A':0, 'C':1, 'G': 2, 'T':3} 
    for seq in dnamatrix: 
    for i in xrange(len(dnamatrix[0])):    
     ProfileMatrix[indices[seq[i]]][i] += 1 
    ProbMatrix = [[float(x)/sum(zip(*ProfileMatrix)[0]) for x in y] for y in ProfileMatrix] 
    return ProbMatrix 
def ProfileRandomGenerator(profile, dna, k, i): 
    indices = {'A':0, 'C':1, 'G': 2, 'T':3} 
    score_list = [] 
    for x in xrange(len(dna[i]) - k + 1): 
     probability = 1 
     window = dna[i][x : k + x] 
    for y in xrange(k): 
     probability *= profile[indices[window[y]]][y] 
    score_list.append(probability) 
    rnd = uniform(0, sum(score_list)) 
    current = 0 
    for z, bias in enumerate(score_list): 
     current += bias 
     if rnd <= current: 
      return dna[i][z : k + z] 
def score(motifs): 
    ProfileMatrix = [[0 for x in xrange(len(motifs[0]))] for x in xrange(4)] 
    indices = {'A':0, 'C':1, 'G': 2, 'T':3} 
    for seq in motifs: 
     for i in xrange(len(motifs[0])):    
      ProfileMatrix[indices[seq[i]]][i] += 1 
    score = len(motifs)*len(motifs[0]) - sum([max(x) for x in zip(*ProfileMatrix)]) 
    return score 
from random import randint, uniform  
def GibbsSampler(k, t, N): 
    dna = ['CGCCCCTCTCGGGGGTGTTCAGTAACCGGCCA', 
    'GGGCGAGGTATGTGTAAGTGCCAAGGTGCCAG', 
    'TAGTACCGAGACCGAAAGAAGTATACAGGCGT', 
    'TAGATCAAGTTTCAGGTGCACGTCGGTGAACC', 
    'AATCCACCAGCTCCACGTGCAATGTTGGCCTA'] 
    Motifs = [] 
    for i in [randint(0, len(dna[0])-k) for x in range(len(dna))]: 
     j = 0 
     kmer = dna[j][i : k+i] 
     j += 1 
     Motifs.append(kmer) 
    BestMotifs = [] 
    s_best = float('inf') 
    for i in xrange(N): 
     x = randint(0, t-1) 
    Motifs.pop(x) 
    profile = BuildProfileMatrix(Motifs) 
    Motif = ProfileRandomGenerator(profile, dna, k, x) 
    Motifs.append(Motif) 
    s_motifs = score(Motifs) 
    if s_motifs < s_best: 
     s_best = s_motifs 
     BestMotifs = Motifs 
return [s_best, BestMotifs] 

k, t, N =8, 5, 100    
best_motifs = [float('inf'), None] 

# Repeat the Gibbs sampler search 20 times. 
for repeat in xrange(20): 
    current_motifs = GibbsSampler(k, t, N) 
    if current_motifs[0] < best_motifs[0]: 
     best_motifs = current_motifs 
# Print and save the answer. 
print '\n'.join(best_motifs[1])    

不幸的是，我的代码从未给出与解决的例子相同的输出。此外，在尝试调试代码时，我发现我得到了定义图案之间不匹配的奇怪分数。但是，当我试图单独运行评分函数时，它完美运行。

每次我运行该脚本时，输出的变化，但无论如何，在这里是存在于代码输出的输入中的一个的示例：

的我的代码

实施例输出

TATGTGTA 
TATGTGTA 
TATGTGTA 
GGTGTTCA 
TATACAGG 

你能帮我调试这段代码吗？我花了整整一天的时间试图找出它有什么问题，尽管我知道这可能是我犯的一个愚蠢的错误，但我的目光未能抓住它。

谢谢大家！

Answer 1

最后，我发现我的代码出了什么问题！这是在第54行：

Motifs.append(Motif) 

后随机删除的主题之一，其次是建立一个轮廓这些图案，然后随机选择基于此配置文件一个新的主题，我应该补充在选定的主题删除前的相同位置不附加到主题列表的末尾。现在

，正确的代码是：

Motifs.insert(x, Motif) 

新的代码工作正常。