Ricerca di motivi con campionatore Gibbs

Question

Sono un principiante sia in programmazione che in bioinformatica. Quindi, apprezzerei la tua comprensione. Ho provato a sviluppare uno script python per la ricerca di motivi usando il campionamento di Gibbs come spiegato nella classe Coursera, "Trovare messaggi nascosti nel DNA". Lo pseudo codice fornito nel corso è:

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 

Descrizione del problema:

codice challenge: Implementare GIBBSSAMPLER.

Ingresso: Numeri interi k, t e N, seguiti da una raccolta di stringhe Dna. Output: le stringhe BestMotif risultanti dall'esecuzione di GIBBSSAMPLER (Dna, k, t, N) con 20 avviamenti casuali. Ricordati di usare pseudocounts!

Esempio Ingresso:

8 5 100 
CGCCCCTCTCGGGGGTGTTCAGTAACCGGCCA 
GGGCGAGGTATGTGTAAGTGCCAAGGTGCCAG 
TAGTACCGAGACCGAAAGAAGTATACAGGCGT 
TAGATCAAGTTTCAGGTGCACGTCGGTGAACC 
AATCCACCAGCTCCACGTGCAATGTTGGCCTA 

Esempio di output:

TCTCGGGG 
CCAAGGTG 
TACAGGCG 
TTCAGGTG 
TCCACGTG 

ho seguito il pseudocodice al meglio delle mie conoscenze. Ecco il mio codice:

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

Sfortunatamente, il mio codice non fornisce mai la stessa uscita dell'esempio risolto. Inoltre, mentre cercavo di eseguire il debug del codice, ho scoperto che ricevo strani punteggi che definiscono le discrepanze tra i motivi. Tuttavia, quando ho provato a eseguire separatamente la funzione score, ha funzionato perfettamente.

Ogni volta che viene eseguita la scrittura, le modifiche di uscita, ma comunque qui è un esempio di una delle uscite per l'ingresso presenti nel codice:

uscita

Esempio del mio codice

TATGTGTA 
TATGTGTA 
TATGTGTA 
GGTGTTCA 
TATACAGG 

Potresti per favore aiutarmi a eseguire il debug di questo codice? !! Ho passato l'intera giornata a cercare di scoprire cosa c'è che non va, anche se so che potrebbe essere un errore stupido che ho fatto, ma il mio occhio non è riuscito a prenderlo.

Grazie a tutti !!

Answer 1

Infine, ho scoperto cosa c'era di sbagliato nel mio codice! Fu in linea 54:

Motifs.append(Motif) 

Dopo aver rimosso in modo casuale uno dei motivi, seguita da costruire un profilo su questi motivi quindi selezionando casualmente un nuovo motivo basato su questo profilo, avrei aggiunto il motivo selezionato nel stessa posizione prima della rimozione NON aggiunta alla fine della lista dei motivi.

Ora, il codice corretto è:

Motifs.insert(x, Motif) 

Il nuovo codice ha funzionato come previsto.