Tìm kiếm motif với lấy mẫu Gibbs

Question

Tôi là người mới bắt đầu trong cả lập trình và tin sinh học. Vì vậy, tôi sẽ đánh giá cao sự hiểu biết của bạn. Tôi đã cố gắng phát triển một kịch bản python cho tìm kiếm motif bằng cách sử dụng lấy mẫu Gibbs như được giải thích trong lớp Coursera, "Tìm kiếm các thông điệp ẩn trong DNA". Các giả cung cấp trong khóa học là:

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 

Mô tả sự cố:

MÃ CHALLENGE: Thực hiện GIBBSSAMPLER.

Nhập: Số nguyên k, t và N, theo sau là tập hợp các chuỗi Dna. Đầu ra: Các chuỗi BestMotifs phát sinh từ việc chạy GIBBSSAMPLER (Dna, k, t, N) với 20 bắt đầu ngẫu nhiên. Hãy nhớ sử dụng các giả tưởng!

Sample Input:

8 5 100 
CGCCCCTCTCGGGGGTGTTCAGTAACCGGCCA 
GGGCGAGGTATGTGTAAGTGCCAAGGTGCCAG 
TAGTACCGAGACCGAAAGAAGTATACAGGCGT 
TAGATCAAGTTTCAGGTGCACGTCGGTGAACC 
AATCCACCAGCTCCACGTGCAATGTTGGCCTA 

Sample Output:

TCTCGGGG 
CCAAGGTG 
TACAGGCG 
TTCAGGTG 
TCCACGTG 

Tôi làm theo các giả đến theo sự hiểu biết của tôi. Đây là mã của tôi:

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

Thật không may, mã của tôi không bao giờ đưa ra kết quả tương tự như ví dụ đã được giải quyết. Bên cạnh đó, trong khi cố gắng để gỡ lỗi mã tôi thấy rằng tôi nhận được điểm số kỳ lạ xác định sự không phù hợp giữa các họa tiết. Tuy nhiên, khi tôi cố gắng chạy hàm số điểm riêng biệt, nó hoạt động hoàn hảo.

Mỗi lần tôi chạy kịch bản, những thay đổi đầu ra, nhưng dù sao đây là một ví dụ về một trong những kết quả đầu ra cho các đầu vào hiện nay trong các mã:

Ví dụ đầu ra của mã của tôi

TATGTGTA 
TATGTGTA 
TATGTGTA 
GGTGTTCA 
TATACAGG 

Bạn có thể giúp tôi gỡ lỗi mã này không? !! Tôi đã dành cả ngày cố gắng tìm ra những gì sai với nó mặc dù tôi biết nó có thể là một số sai lầm ngớ ngẩn mà tôi đã thực hiện, nhưng mắt tôi không bắt được nó.

Cảm ơn tất cả !!

Answer 1

Cuối cùng, tôi đã tìm ra điều gì sai trong mã của mình! Đó là vào dòng 54:

Motifs.append(Motif) 

Sau khi loại bỏ một cách ngẫu nhiên một trong những họa tiết, tiếp theo là xây dựng một hồ sơ cá nhân ra khỏi những mô típ sau đó chọn ngẫu nhiên một motif mới dựa trên hồ sơ này, tôi nên đã thêm các motif chọn trong cùng một vị trí trước khi xóa KHÔNG được nối vào cuối danh sách motif.

Bây giờ, mã đúng là:

Motifs.insert(x, Motif) 

Mã mới làm việc như mong đợi.