title: Motif search with Gibbs sampler
body: <p>I am a beginner in both programming and bioinformatics. So, I would appreciate your understanding. I tried to develop a python script for motif search using Gibbs sampling as explained in Coursera class, "Finding Hidden Messages in DNA". The pseudocode provided in the course is:</p>
  
  <pre><code>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) &lt; Score(BestMotifs)
              BestMotifs ← Motifs
      return BestMotifs
  </code></pre>
  
  <p>Problem description:</p>
  
  <p>CODE CHALLENGE: Implement GIBBSSAMPLER.</p>
  
  <p><strong>Input:</strong> Integers k, t, and N, followed by a collection of strings Dna.
       Output: The strings BestMotifs resulting from running GIBBSSAMPLER(Dna, k, t, N) with
       20 random starts. Remember to use pseudocounts!</p>
  
  <p><strong>Sample Input</strong>:      </p>
  
  <pre><code> 8 5 100
   CGCCCCTCTCGGGGGTGTTCAGTAACCGGCCA
   GGGCGAGGTATGTGTAAGTGCCAAGGTGCCAG
   TAGTACCGAGACCGAAAGAAGTATACAGGCGT
   TAGATCAAGTTTCAGGTGCACGTCGGTGAACC
   AATCCACCAGCTCCACGTGCAATGTTGGCCTA
  </code></pre>
  
  <p><strong>Sample Output</strong>:</p>
  
  <pre><code> TCTCGGGG
   CCAAGGTG
   TACAGGCG
   TTCAGGTG
   TCCACGTG
  </code></pre>
  
  <p>I followed the pseudocode to the best of my knowledge. Here is my code:</p>
  
  <pre><code>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 &lt;= 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 &lt; 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] &lt; best_motifs[0]:
          best_motifs = current_motifs
  # Print and save the answer.
  print '\n'.join(best_motifs[1])            
  </code></pre>
  
  <p>Unfortunately, my code never gives the same output as the solved example. Besides, while trying to debug the code I found that I get weird scores that define the mismatches between motifs. However, when I tried to run the score function separately, it worked perfectly.</p>
  
  <p>Each time I run the script, the output changes, but anyway here is an example of one of the outputs for the input present in the code:</p>
  
  <p><strong>Example output of my code</strong></p>
  
  <pre><code>TATGTGTA
  TATGTGTA
  TATGTGTA
  GGTGTTCA
  TATACAGG
  </code></pre>
  
  <p>Could you please help me debug this code?!! I spent the whole day trying to find out what's wrong with it although I know it might be some silly mistake I made, but my eye failed to catch it.</p>
  
  <p>Thank you all!!</p>