I start both in programming and in bioinformatics. So, I would appreciate your understanding. I tried to develop a python script to search for motifs using a Gibbs selection, as described in the Coursera class, "Search for hidden messages in DNA." The pseudocode indicated in the course is as follows:
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
Description of the problem:
CODE CHALLENGE: Implementing the GIBBSSAMPLER.
Input: Integers k, t, and N, followed by a collection of Dna strings. Result: BestMotifs strings resulting from running GIBBSSAMPLER (Dna, k, t, N) with 20 random starts. Remember to use pseudo packages!
Input Example :
8 5 100 CGCCCCTCTCGGGGGTGTTCAGTAACCGGCCA GGGCGAGGTATGTGTAAGTGCCAAGGTGCCAG TAGTACCGAGACCGAAAGAAGTATACAGGCGT TAGATCAAGTTTCAGGTGCACGTCGGTGAACC AATCCACCAGCTCCACGTGCAATGTTGGCCTA
Output result :
TCTCGGGG CCAAGGTG TACAGGCG TTCAGGTG TCCACGTG
I followed the pseudo code as far as I know. Here is my 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 <= 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])
Unfortunately, my code never gives the same result as the solved example. In addition, when I tried to debug the code, I found that I was getting strange ratings that determined inconsistencies between the motives. However, when I tried to run the evaluation function separately, it worked fine.
Each time I run the script, the output changes, but in any case, here is an example of one of the outputs for input present in the code:
Sample output of my code
TATGTGTA TATGTGTA TATGTGTA GGTGTTCA TATACAGG
Could you help me debug this code? !! I spent the whole day trying to figure out what was wrong with him, although I know that it could be some kind of stupid mistake that I made, but my eye could not catch it.
Thanks everyone!