使用 Gibbs 采样器进行 Motif 搜索

Question

我是编程和生物信息学的初学者。所以，我会很感激你的理解。我尝试使用 Gibbs 采样开发一个用于主题搜​​索的 python 脚本，如 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

问题描述：

代码挑战：实施 GIBBSSAMPLER。

输入： 整数 k、t 和 N，后跟字符串 Dna 的集合。 输出：运行 GIBBSSAMPLER(Dna, k, t, N) 产生的字符串 BestMotifs 20个随机开始。记得使用伪计数！

示例输入：

 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])            

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

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

我的代码的示例输出

TATGTGTA
TATGTGTA
TATGTGTA
GGTGTTCA
TATACAGG

你能帮我调试一下这段代码吗？！！我花了一整天的时间试图找出它出了什么问题，虽然我知道这可能是我犯的一些愚蠢的错误，但我的眼睛没能抓住它。

谢谢大家！！

Answer 1

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

Motifs.append(Motif)

在随机删除一个主题后，然后从这些主题中构建一个配置文件，然后根据该配置文件随机选择一个新主题，我应该在删除之前将选定的主题添加到相同的位置，而不是附加到末尾主题列表。

现在，正确的代码是：

Motifs.insert(x, Motif)

新代码按预期工作。