python - 使用 Gibbs 采样器进行基序搜索

Question

我是编程和生物信息学的初学者。因此，非常感谢您的理解。我尝试使用 Gibbs 采样开发一个用于主题搜​​索的 python 脚本，如 Coursera 类(class)“Finding Hidden Messages in DNA”中所述。类(class)中提供的伪代码是:

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

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

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

我的代码的示例输出

TATGTGTA
TATGTGTA
TATGTGTA
GGTGTTCA
TATACAGG

你能帮我调试这段代码吗？!!我花了一整天时间试图找出问题所在，虽然我知道这可能是我犯的一些愚蠢的错误，但我的眼睛却没有注意到。

谢谢大家!!

Answer 1

终于，我发现我的代码出了什么问题!它在第 54 行:

Motifs.append(Motif)

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

现在，正确的代码是:

Motifs.insert(x, Motif)

新代码按预期工作。