<commit_before>#!/usr/bin/python
import argparse
import itertools
import random
import copy
from collections import Counter

class HiddenMessages:
    def __init__(self):
        self.author = '<NAME>'
        self.substitution = {
            'A': 'T',
            'T': 'A',
            'G': 'C',
            'C': 'G'
        }
        self.lex_order = {
            'A': 0,
            'C': 1,
            'G': 2,
            'T': 3
        }
        self.reverse_lex_order = {
            '0': 'A',
            '1': 'C',
            '2': 'G',
            '3': 'T'
        }

    def gibbs_sampler(self, k, t, N, dna):
        motifs = []
        for dnai in dna:
            i = random.randint(0, len(dnai) - k)
            motifs.append(dnai[i:i+k])

        best_motifs = copy.deepcopy(motifs) # Deep copy
        for n in xrange(N):
            i = random.randint(0, len(motifs) - 1)
            profiles = hm.create_laplace_rule_profile([motifs[j] for j in xrange(len(motifs)) if j != i])
            # print [motifs[j] for j in xrange(len(motifs)) if j != i]
            # print profiles

            # Calculate probabilty of all k-mers in dnai
            all_kmers_in_dnai = [dna[i][a:a+k] for a in xrange(len(dna[i]) - k + 1)]
            prob_dist = [reduce(lambda x,y:x*y, [profiles[kmer[nuci]][nuci] for nuci in xrange(len(kmer))]) for kmer in all_kmers_in_dnai]
            # print prob_dist
            new_motif_i = hm.random_biased_die(prob_dist)
            motifs[i] = all_kmers_in_dnai[new_motif_i - 1]  # Because random biased die always starts at 1

            # print best_motifs
            # print motifs
            # print 'old: %s new: %s'%(hm.score_motifs(best_motifs), hm.score_motifs(motifs))
            if hm.score_motifs(motifs) < hm.score_motifs(best_motifs):
                best_motifs = copy.deepcopy(motifs) # Deep copy

        print hm.score_motifs(best_motifs)
        return best_motifs

    def random_biased_die(self, prob_dist):
        prob_dist_to_one = [i/sum(prob_dist) for i in prob_dist]
        cumulative_prob_dist_to_one = [sum(prob_dist_to_one[:i+1]) for i in range(len(prob_dist_to_one))]

        num = random.random()
        for i in xrange(len(prob_dist)):
            if num < cumulative_prob_dist_to_one[i]:
                return i + 1

        return len(prob_dist) + 1

    def randomized_motif_search(self, k , t, dna):
        def inner(k, t, dna):
            motifs = []
            for dnai in dna:
                i = random.randint(0, len(dnai) - k)
                motifs.append(dnai[i:i+k])

            while True:
                profiles = hm.create_laplace_rule_profile(motifs)
                motif_prime = [hm.profile_most_probable_kmer(dnai, k, profiles['A'], profiles['C'], profiles['G'], profiles['T']) for dnai in dna]

                if hm.score_motifs(motif_prime) < hm.score_motifs(motifs):
                    motifs = motif_prime
                else:
                    return motifs

        best_motifs = []
        for i in xrange(1000):
            # best_motifs = inner(k, t, dna)
            if not best_motifs:
                best_motifs = inner(k, t, dna)
            else:
                temp_motifs = inner(k, t, dna)
                if hm.score_motifs(temp_motifs) < hm.score_motifs(best_motifs):
                    best_motifs = temp_motifs

        print hm.score_motifs(best_motifs)
        return best_motifs

    def greedy_motif_search_w_pseudocounts(self, dna, k, t):
        best_motifs = [dnai[:k] for dnai in dna]
        score_best_motifs = hm.score_motifs(best_motifs)

        for i in xrange(len(dna[0]) - k + 1):
            motif = []
            motif.append(dna[0][i:i+k])
            for j in xrange(1, t):
                profiles = hm.create_laplace_rule_profile(motif)
                motifj = hm.profile_most_probable_kmer(dna[j], k, profiles['A'], profiles['C'], profiles['G'], profiles['T'])
                motif.append(motifj)

            motifs = motif

            if hm.score_motifs(motifs) < score_best_motifs:
                best_motifs = motifs
                score_best_motifs = hm.score_motifs(motifs)

        return best_motifs

    def create_laplace_rule_profile(self, motifs):
        k = len(motifs[0])
        profile = {
            'A': [0]*k,
            'C': [0]*k,
            'G': [0]*k,
            'T': [0]*k,
        }

        for nuc in profile:
            for i in xrange(k):
                # profile[nuc][i] = sum([2 if motif[i] == nuc else 1 for motif in motifs])*1.0/len(motifs)
                profile[nuc][i] = (sum([1 if motif[i] == nuc else 0 for motif in motifs])+1)*1.0/len(motifs)

        return profile

    def greedy_motif_search(self, dna, k, t):
        best_motifs = [dnai[:k] for dnai in dna]
        score_best_motifs = hm.score_motifs(best_motifs)

        for i in xrange(len(dna[0]) - k + 1):
            motif = []
            motif.append(dna[0][i:i+k])
            for j in xrange(1, t):
                profiles = hm.create_profile(motif)
                motifj = hm.profile_most_probable_kmer(dna[j], k, profiles['A'], profiles['C'], profiles['G'], profiles['T'])
                motif.append(motifj)

            motifs = motif

            if hm.score_motifs(motifs) < score_best_motifs:
                best_motifs = motifs
                score_best_motifs = hm.score_motifs(motifs)

        return best_motifs

    def create_profile(self, motifs):
        k = len(motifs[0])
        profile = {
            'A': [0]*k,
            'C': [0]*k,
            'G': [0]*k,
            'T': [0]*k,
        }

        for nuc in profile:
            for i in xrange(k):
                profile[nuc][i] = sum([1 if motif[i] == nuc else 0 for motif in motifs])*1.0/len(motifs)

        return profile

    def score_motifs(self, motifs):
        transposed = [''.join([motif[i] for motif in motifs]) for i in xrange(len(motifs[0]))]

        return sum([len(i) - Counter(i).most_common(1)[0][1] for i in transposed])


    def profile_most_probable_kmer(self, dna, k, pA, pC, pG, pT):
        profiles = {
            'A': pA,
            'C': pC,
            'G': pG,
            'T': pT
        }
        highest_prob = -1 
        highest_prob_kmer = ''

        for i in xrange(len(dna) - k + 1):
            current_kmer = dna[i:i+k]
            current_kmer_prob = reduce(lambda a, b: a*b, [profiles[current_kmer[nuc_i]][nuc_i] for nuc_i in xrange(len(current_kmer))])

            if current_kmer_prob > highest_prob:
                highest_prob = current_kmer_prob
                highest_prob_kmer = current_kmer

        return highest_prob_kmer

    def all_median_strings(self, dna, k):
        distance = 9999999999
        median = []

        for pattern in itertools.product('ACTG', repeat=k):
            pattern = ''.join(pattern)
            d_pattern_dna = sum([min([self.hamming_distance(pattern, dna_i[j:j+k]) for j in range(len(dna_i) - k + 1)]) for dna_i in dna])

            if distance == d_pattern_dna:
                median.append(pattern)
            elif distance > d_pattern_dna:
                distance = d_pattern_dna
                median = [pattern]

        return median

    def median_string(self, dna, k):
        distance = 9999999999
        median = ''

        for pattern in itertools.product('ACTG', repeat=k):
            pattern = ''.join(pattern)
            d_pattern_dna = sum([min([self.hamming_distance(pattern, dna_i[j:j+k]) for j in range(len(dna_i) - k + 1)]) for dna_i in dna])

            if distance > d_pattern_dna:
                distance = d_pattern_dna
                median = pattern

        return median

    def motif_enumeration(self, dna_set, k, d):
        patterns = set()
        raw_patterns = set()
        for dna in dna_set:
            for raw_pattern in [dna[i:i+k] for i in range(len(dna) - k + 1)]:
                raw_patterns.add(raw_pattern)

        for pattern in raw_patterns:
            for pattern_prime in self.neighbors(pattern, d):
                appears = True
                for dna in dna_set:
                    if self.approx_pattern_count(pattern_prime, dna, d) == 0:
                        appears = False
                        break

                if appears:
                    patterns.add(pattern_prime)

        return patterns

    # From: http://stackoverflow.com/a/19941659
    def neighbors(self, word, hamming_distance, charset='ATCG'):
        for indices in itertools.combinations(range(len(word)), hamming_distance):
            for replacements in itertools.product(charset, repeat=hamming_distance):
                mutation = list(word)
                for index, replacement in zip(indices, replacements):
                    mutation[index] = replacement
                yield "".join(mutation)

    # This is the brute force approach
    def freq_words_w_mismatches_and_reverse_comp(self, text, k, d):
        freq = [0] * (4**k)

        # Calculate freq
        for i in xrange(4**k):
            pattern = self.number_to_pattern(i, k)
            freq[i] = self.approx_pattern_count(pattern, text, d) + self.approx_pattern_count(self.reverse_complement(pattern), text, d)

        # Return only max words
        maxFreq = max(freq)
        maxWords = []
        for i in xrange(4**k):
            if freq[i] == maxFreq:
                maxWords.append(self.number_to_pattern(i, k))

        return maxWords

    # This is the brute force approach
    def freq_words_w_mismatches(self, text, k , d):
        freq = [0] * (4**k)

        # Calculate freq
        for i in xrange(4**k):
            freq[i] = self.approx_pattern_count(self.number_to_pattern(i, k), text, d)

        # Return only max words
        maxFreq = max(freq)
        maxWords = []
        for i in xrange(4**k):
            if freq[i] == maxFreq:
                maxWords.append(self.number_to_pattern(i, k))

        return maxWords

    def approx_pattern_count(self, pattern, text, d):
        return len(self.approx_pattern_matching(pattern, text, d))

    def approx_pattern_matching(self, pattern, text, d):
        len_pattern = len(pattern)
        pos = []
        for i in xrange(len(text) - len_pattern + 1):
            if self.hamming_distance(text[i:i+len_pattern], pattern) <= d:
                pos.append(i)

        return pos

    def hamming_distance(self, p, q):
        hamming_distance = 0
        for i in xrange(len(p)):
            if p[i] != q[i]:
                hamming_distance += 1

        return hamming_distance

    def min_skew(self, genome):
        min_skew_value = 0
        min_skew_i = []
        skew_array = [0]

        for i in xrange(len(genome)):
            base = genome[i]
            if base == 'C':
                skew_array.append(skew_array[-1] - 1)
            elif base == 'G':
                skew_array.append(skew_array[-1] + 1)
            else:
                skew_array.append(skew_array[-1])

            if skew_array[-1] < min_skew_value:
                min_skew_value = skew_array[-1]
                min_skew_i = [i + 1]
            elif skew_array[-1] == min_skew_value:
                min_skew_i.append(i + 1)

        return min_skew_i 

    def skew_genome(self, genome):
        skew_array = [0]

        for base in genome:
            if base == 'C':
                skew_array.append(skew_array[-1] - 1)
            elif base == 'G':
                skew_array.append(skew_array[-1] + 1)
            else:
                skew_array.append(skew_array[-1])

        return skew_array

    def pattern_to_number2(self, pattern):
        if not pattern:
            return 0

        symbol = pattern[-1]
        prefix = pattern[:-1]
        return 4 * self.pattern_to_number2(prefix) + self.lex_order[symbol]

    def computing_frequencies(self, text, k):
        freq = [0] * (4**k)
        for i in xrange(len(text) - k + 1):
            pattern = text[i:i+k]
            j = self.pattern_to_number(pattern)
            #print '%s -> %s'%(pattern, j)
            freq[j] += 1

        return freq

    def number_to_pattern(self, index, k):
        # 11, 2 = GT
        def recurse(num):
            if num < 4:
                return str(num)
            else:
                return str(recurse(num/4)) + str(num%4)

        base4_in_dna = [self.reverse_lex_order[i] for i in recurse(index)]
        return ''.join(['A']*(k-len(base4_in_dna)) + base4_in_dna)
        

    def pattern_to_number(self, pattern):
        # This is essentially base 4
        # GT = 11
        number = 0
        for i in xrange(len(pattern)):
            number += 4**(len(pattern) - 1 - i) * self.lex_order[pattern[i]]

        return number

    def clump_finder2(self, genome, k, L, t):
        # Generate freq array 
        freq_array = {}

        for i in range(len(genome) - k + 1):
            if genome[i:i+k] not in freq_array:
                freq_array[genome[i:i+k]] = [i]
            else:
                freq_array[genome[i:i+k]].append(i)

        # Calculate Lt-clumps
        tWords = []
        for word in freq_array:
            if len(freq_array[word]) < t:
                continue
            else:
                for i in range(len(freq_array[word]) - t + 1):
                    if freq_array[word][i+t-1] + k - freq_array[word][i] <= L:
                        tWords.append(word)
                        break

        return tWords
        

    def clump_finder(self, genome, k, L, t):
        def frequent_words4(text, k, t):
            # Do everything in one pass
            wordCounts = {}
            tWords = {};
        
            for i in range(len(text) - k + 1):
                currWord = text[i:i+k]
                
                # Update freq count of word in currWords dict
                if currWord in wordCounts:
                    wordCounts[currWord] += 1
                else:
                    wordCounts[currWord] = 1
        
                # See whether this is a new tWord entry
                if wordCounts[currWord] >= t:
                    tWords[currWord] = True;
        
            return tWords.keys()

        kLt_patterns = []
        for i in xrange(len(genome) - L):
            window_genome = genome[i:i+L]
            new_patterns = frequent_words4(window_genome, k, t)
            kLt_patterns += new_patterns

        return set(kLt_patterns)

    def pattern_matching(self, pattern, genome):
        patternLength = len(pattern)
        positions = []
    
        for i in range(len(genome)):
            if genome[i:i+patternLength] == pattern:
                positions.append(i)
    
        return ' '.join([str(i) for i in positions])

    def reverse_complement(self, pattern):
        return ''.join([self.substitution[i] for i in pattern][::-1])

    def frequent_words(self, text, k):
        # Generate all words of len k from text
        words = {}
        for i in range(len(text) - k + 1):
            words[text[i:i+k]] = 0
    
        # Count freq of words
        # and add entry into wordsByFreq as well
        wordsByFreq = {}
        for word in words.keys():
            words[word] = self.pattern_count(text, word)
    
            if words[word] not in wordsByFreq:
                wordsByFreq[words[word]] = [word]
            else:
                wordsByFreq[words[word]].append(word)
    
        # Return maximum
        return ' '.join(wordsByFreq[max(wordsByFreq.keys())])
    
    def frequent_words2(self, text, k):
        # Do everything in one pass
        wordCounts = {}
    
        for i in range(len(text) - k + 1):
            currWord = text[i:i+k]
            if currWord in wordCounts:
                wordCounts[currWord] += 1
            else:
                wordCounts[currWord] = 1
    
        maxFreq = max(wordCounts.values())
        return ' '.join([i for i in wordCounts.keys() if wordCounts[i] == maxFreq])
    
    def frequent_words3(self, text, k):
        # Do everything in one pass
        wordCounts = {}
        currMax = 0;
        maxWords = [];
    
        for i in range(len(text) - k + 1):
            currWord = text[i:i+k]
    
            # Update freq count of word in currWords dict
            if currWord in wordCounts:
                wordCounts[currWord] += 1
            else:
                wordCounts[currWord] = 1
    
            # See whether this is a new maxWord entry
            if wordCounts[currWord] == currMax:
                maxWords.append(currWord)
            elif wordCounts[currWord] > currMax:
                currMax = wordCounts[currWord]
                maxWords = [currWord]
    
        return ' '.join(maxWords)

    def pattern_count(self, text, pattern):
        patternLength = len(pattern)
        count = 0
    
        for i in range(len(text)):
            if text[i:i+patternLength] == pattern:
                count += 1
    
        return count

if __name__ == "__main__":
    hm = HiddenMessages()

    parser = argparse.ArgumentParser(description='Coursera Hidden Messages - Jhalley')
    parser.add_argument('--pattern_count', help='Count the number of occurances of pattern in text')
    parser.add_argument('--freq_words', help='Find most frequent k-mers in text')
    parser.add_argument('--freq_words2', help='Find most frequent k-mers in text - faster')
    parser.add_argument('--freq_words3', help='Find most frequent k-mers in text - fastest!')
    parser.add_argument('--reverse_complement', help='Find the reverse complement of a DNA string')
    parser.add_argument('--pattern_matching', help='Find all occurrences of a pattern in a string')
    parser.add_argument('--clump_finder', help='Find patterns forming clumps in a string')
    parser.add_argument('--clump_finder2', help='Find patterns forming clumps in a string using freq array')
    parser.add_argument('--computing_frequencies', help='Generate frequency array')
    parser.add_argument('--pattern_to_number2', help='Pattern to number algorithm')
    parser.add_argument('--number_to_pattern', help='My number to pattern algorithm')
    parser.add_argument('--skew_genome', help='Get skew diagram of genome')
    parser.add_argument('--min_skew', help='Find a position in a genome where the skew diagram attains a minimum')
    parser.add_argument('--hamming_distance', help='Calculate hamming distance between two strings of equal length')
    parser.add_argument('--approx_pattern_matching', help='Find all approximate occurrences of a pattern in a string')
    parser.add_argument('--approx_pattern_count', help='Get count of all approximate occurrences of a pattern in a string')
    parser.add_argument('--freq_words_w_mismatches', help='Find the most frequent k-mers with mismatches in a string')
    parser.add_argument('--freq_words_w_mismatches_and_reverse_comp', help='Find the most frequent k-mers (with mismatches and reverse complements) in a string')
    parser.add_argument('--motif_enumeration', help='A brute force attempt at motif finding')
    parser.add_argument('--median_string', help='Find a median string')
    parser.add_argument('--profile_most_probable_kmer', help='Find a Profile-most probable k-mer in a string')
    parser.add_argument('--greedy_motif_search', help='Greedy motif search')
    parser.add_argument('--greedy_motif_search_w_pseudocounts', help='Greedy motif search with pseudocounts')
    parser.add_argument('--randomized_motif_search', help='Randomized motif search')
    parser.add_argument('--gibbs_sampler', help='Gibbs sampler')
    args = parser.parse_args()

     
    # Parse from file
    if args.pattern_count:
        with open(args.pattern_count) as f:
            lines = f.readlines()
        print hm.pattern_count(lines[0].strip(), lines[1].strip())
    elif args.freq_words:
        with open(args.freq_words) as f:
            lines = f.readlines()
        print hm.frequent_words(lines[0].strip(), int(lines[1].strip()))
    elif args.freq_words2:
        with open(args.freq_words2) as f:
            lines = f.readlines()
        print hm.frequent_words2(lines[0].strip(), int(lines[1].strip()))
    elif args.freq_words3:
        with open(args.freq_words3) as f:
            lines = f.readlines()
        print hm.frequent_words3(lines[0].strip(), int(lines[1].strip()))
    elif args.reverse_complement:
        with open(args.reverse_complement) as f:
            lines = f.readlines()
        print hm.reverse_complement(lines[0].strip())
    elif args.pattern_matching:
        with open(args.pattern_matching) as f:
            lines = f.readlines()
        print hm.pattern_matching(lines[0].strip(), lines[1].strip())
    elif args.clump_finder:
        with open(args.clump_finder) as f:
            lines = f.readlines()
        klt = [int(i) for i in lines[1].strip().split()]
        print ' '.join(hm.clump_finder(lines[0].strip(), klt[0], klt[1], klt[2]))
    elif args.clump_finder2:
        with open(args.clump_finder2) as f:
            lines = f.readlines()
        klt = [int(i) for i in lines[1].strip().split()]
        print ' '.join(hm.clump_finder2(lines[0].strip(), klt[0], klt[1], klt[2]))
    elif args.computing_frequencies:
        with open(args.computing_frequencies) as f:
            lines = f.readlines()
        print ' '.join([str(i) for i in hm.computing_frequencies(lines[0].strip(), int(lines[1].strip()))])
    elif args.pattern_to_number2:
        with open(args.pattern_to_number2) as f:
            lines = f.readlines()
        print hm.pattern_to_number2(lines[0].strip())
    elif args.number_to_pattern:
        with open(args.number_to_pattern) as f:
            lines = f.readlines()
        print hm.number_to_pattern(int(lines[0].strip()), int(lines[1].strip()))
    elif args.skew_genome:
        with open(args.skew_genome) as f:
            lines = f.readlines()
        print hm.skew_genome(lines[0].strip())
    elif args.min_skew:
        with open(args.min_skew) as f:
            lines = f.readlines()
        print ' '.join([str(i) for i in hm.min_skew(lines[0].strip())])
    elif args.hamming_distance:
        with open(args.hamming_distance) as f:
            lines = f.readlines()
        print hm.hamming_distance(lines[0].strip(), lines[1].strip())
    elif args.approx_pattern_matching:
        with open(args.approx_pattern_matching) as f:
            lines = f.readlines()
        print ' '.join([str(i) for i in hm.approx_pattern_matching(lines[0].strip(), lines[1].strip(), int(lines[2].strip()))])
    elif args.approx_pattern_count:
        with open(args.approx_pattern_count) as f:
            lines = f.readlines()
        print hm.approx_pattern_count(lines[0].strip(), lines[1].strip(), int(lines[2].strip()))
    elif args.freq_words_w_mismatches:
        with open(args.freq_words_w_mismatches) as f:
            lines = f.readlines()
        print ' '.join(hm.freq_words_w_mismatches(lines[0].strip(), int(lines[1].strip().split()[0]), int(lines[1].strip().split()[1])))
    elif args.freq_words_w_mismatches_and_reverse_comp:
        with open(args.freq_words_w_mismatches_and_reverse_comp) as f:
            lines = f.readlines()
        print ' '.join(hm.freq_words_w_mismatches_and_reverse_comp(lines[0].strip(), int(lines[1].strip().split()[0]), int(lines[1].strip().split()[1])))
    elif args.motif_enumeration:
        with open(args.motif_enumeration) as f:
            lines = f.readlines()
        print ' '.join(sorted(hm.motif_enumeration([dna.strip() for dna in lines[1:]], int(lines[0].split()[0]), int(lines[0].split()[1]))))
    elif args.median_string:
        with open(args.median_string) as f:
            lines = f.readlines()
        print hm.median_string([dna.strip() for dna in lines[1:]], int(lines[0].strip()))
    elif args.profile_most_probable_kmer:
        with open(args.profile_most_probable_kmer) as f:
            lines = f.readlines()
        print hm.profile_most_probable_kmer(lines[0].strip(), int(lines[1].strip()), [float(i) for i in lines[2].strip().split()], [float(i) for i in lines[3].strip().split()], [float(i) for i in lines[4].strip().split()], [float(i) for i in lines[5].strip().split()])
    elif args.greedy_motif_search:
        with open(args.greedy_motif_search) as f:
            lines = f.readlines()
        print '\n'.join(hm.greedy_motif_search([dna.strip() for dna in lines[1:]], int(lines[0].strip().split()[0]), int(lines[0].strip().split()[1])))
    elif args.greedy_motif_search_w_pseudocounts:
        with open(args.greedy_motif_search_w_pseudocounts) as f:
            lines = f.readlines()
        print '\n'.join(hm.greedy_motif_search_w_pseudocounts([dna.strip() for dna in lines[1:]], int(lines[0].strip().split()[0]), int(lines[0].strip().split()[1])))
    elif args.randomized_motif_search:
        with open(args.randomized_motif_search) as f:
            lines = f.readlines()
        print '\n'.join(hm.randomized_motif_search(int(lines[0].strip().split()[0]), int(lines[0].strip().split()[1]), [dna.strip() for dna in lines[1:]]))
    elif args.gibbs_sampler:
        with open(args.gibbs_sampler) as f:
            lines = f.readlines()
        print '\n'.join(hm.gibbs_sampler(int(lines[0].strip().split()[0]), int(lines[0].strip().split()[1]), int(lines[0].strip().split()[2]), [dna.strip() for dna in lines[1:]]))

    # Test calls
    # print hm.score_motifs(['TCTCGGGG', 'CCAAGGTG', 'TACAGGCG', 'TTCAGGTG', 'TCCACGTG'])
    #print '\n'.join(hm.gibbs_sampler(8, 5, 1000, ['CGCCCCTCTCGGGGGTGTTCAGTAACCGGCCA', 'GGGCGAGGTATGTGTAAGTGCCAAGGTGCCAG', 'TAGTACCGAGACCGAAAGAAGTATACAGGCGT', 'TAGATCAAGTTTCAGGTGCACGTCGGTGAACC', 'AATCCACCAGCTCCACGTGCAATGTTGGCCTA']))
    #print Counter([hm.random_biased_die([0.1, 0.2, 0.3]) for i in xrange(10000)])
    #print hm.randomized_motif_search(8, 5, ['CGCCCCTCTCGGGGGTGTTCAGTAAACGGCCA', 'GGGCGAGGTATGTGTAAGTGCCAAGGTGCCAG', 'TAGTACCGAGACCGAAAGAAGTATACAGGCGT', 'TAGATCAAGTTTCAGGTGCACGTCGGTGAACC', 'AATCCACCAGCTCCACGTGCAATGTTGGCCTA'])
    # print hm.all_median_strings(['CTCGATGAGTAGGAAAGTAGTTTCACTGGGCGAACCACCCCGGCGCTAATCCTAGTGCCC', 'GCAATCCTACCCGAGGCCACATATCAGTAGGAACTAGAACCACCACGGGTGGCTAGTTTC', 'GGTGTTGAACCACGGGGTTAGTTTCATCTATTGTAGGAATCGGCTTCAAATCCTACACAG'], 7)<commit_msg>Modified gibbs sampler to allow us to run it multiple times<commit_after>#!/usr/bin/python
import argparse
import itertools
import random
import copy
from collections import Counter

class HiddenMessages:
    def __init__(self):
        self.author = '<NAME>'
        self.substitution = {
            'A': 'T',
            'T': 'A',
            'G': 'C',
            'C': 'G'
        }
        self.lex_order = {
            'A': 0,
            'C': 1,
            'G': 2,
            'T': 3
        }
        self.reverse_lex_order = {
            '0': 'A',
            '1': 'C',
            '2': 'G',
            '3': 'T'
        }

    def gibbs_sampler(self, k, t, N, dna):
        # This is a modified gibbs sampler where we run the gibbs sampler multiple times to try and get the global best score
        def inner(k, t, N, dna):
            motifs = []
            for dnai in dna:
                i = random.randint(0, len(dnai) - k)
                motifs.append(dnai[i:i+k])

            local_best_motifs = copy.deepcopy(motifs) # Deep copy
            for n in xrange(N):
                i = random.randint(0, len(motifs) - 1)
                profiles = hm.create_laplace_rule_profile([motifs[j] for j in xrange(len(motifs)) if j != i])

                # Calculate probabilty of all k-mers in dnai
                all_kmers_in_dnai = [dna[i][a:a+k] for a in xrange(len(dna[i]) - k + 1)]
                prob_dist = [reduce(lambda x,y:x*y, [profiles[kmer[nuci]][nuci] for nuci in xrange(len(kmer))]) for kmer in all_kmers_in_dnai]
                new_motif_i = hm.random_biased_die(prob_dist)
                motifs[i] = all_kmers_in_dnai[new_motif_i - 1]  # Because random biased die always starts at 1

                if hm.score_motifs(motifs) < hm.score_motifs(local_best_motifs):
                    local_best_motifs = copy.deepcopy(motifs) # Deep copy

            return local_best_motifs

        best_motifs = []
        for i in xrange(1): # Change this number to run gibbs_sampler multiple times to try to get the global best score
            if not best_motifs:
                best_motifs = inner(k, t, N, dna)
            else:
                temp_motifs = inner(k, t, N, dna)
                if hm.score_motifs(temp_motifs) < hm.score_motifs(best_motifs):
                    best_motifs = copy.deepcopy(temp_motifs)

        print hm.score_motifs(best_motifs)
        return best_motifs

    def random_biased_die(self, prob_dist):
        prob_dist_to_one = [i/sum(prob_dist) for i in prob_dist]
        cumulative_prob_dist_to_one = [sum(prob_dist_to_one[:i+1]) for i in range(len(prob_dist_to_one))]

        num = random.random()
        for i in xrange(len(prob_dist)):
            if num < cumulative_prob_dist_to_one[i]:
                return i + 1

        return len(prob_dist) + 1

    def randomized_motif_search(self, k , t, dna):
        def inner(k, t, dna):
            motifs = []
            for dnai in dna:
                i = random.randint(0, len(dnai) - k)
                motifs.append(dnai[i:i+k])

            while True:
                profiles = hm.create_laplace_rule_profile(motifs)
                motif_prime = [hm.profile_most_probable_kmer(dnai, k, profiles['A'], profiles['C'], profiles['G'], profiles['T']) for dnai in dna]

                if hm.score_motifs(motif_prime) < hm.score_motifs(motifs):
                    motifs = motif_prime
                else:
                    return motifs

        best_motifs = []
        for i in xrange(1000):
            # best_motifs = inner(k, t, dna)
            if not best_motifs:
                best_motifs = inner(k, t, dna)
            else:
                temp_motifs = inner(k, t, dna)
                if hm.score_motifs(temp_motifs) < hm.score_motifs(best_motifs):
                    best_motifs = temp_motifs

        print hm.score_motifs(best_motifs)
        return best_motifs

    def greedy_motif_search_w_pseudocounts(self, dna, k, t):
        best_motifs = [dnai[:k] for dnai in dna]
        score_best_motifs = hm.score_motifs(best_motifs)

        for i in xrange(len(dna[0]) - k + 1):
            motif = []
            motif.append(dna[0][i:i+k])
            for j in xrange(1, t):
                profiles = hm.create_laplace_rule_profile(motif)
                motifj = hm.profile_most_probable_kmer(dna[j], k, profiles['A'], profiles['C'], profiles['G'], profiles['T'])
                motif.append(motifj)

            motifs = motif

            if hm.score_motifs(motifs) < score_best_motifs:
                best_motifs = motifs
                score_best_motifs = hm.score_motifs(motifs)

        return best_motifs

    def create_laplace_rule_profile(self, motifs):
        k = len(motifs[0])
        profile = {
            'A': [0]*k,
            'C': [0]*k,
            'G': [0]*k,
            'T': [0]*k,
        }

        for nuc in profile:
            for i in xrange(k):
                # profile[nuc][i] = sum([2 if motif[i] == nuc else 1 for motif in motifs])*1.0/len(motifs)
                profile[nuc][i] = (sum([1 if motif[i] == nuc else 0 for motif in motifs])+1)*1.0/len(motifs)

        return profile

    def greedy_motif_search(self, dna, k, t):
        best_motifs = [dnai[:k] for dnai in dna]
        score_best_motifs = hm.score_motifs(best_motifs)

        for i in xrange(len(dna[0]) - k + 1):
            motif = []
            motif.append(dna[0][i:i+k])
            for j in xrange(1, t):
                profiles = hm.create_profile(motif)
                motifj = hm.profile_most_probable_kmer(dna[j], k, profiles['A'], profiles['C'], profiles['G'], profiles['T'])
                motif.append(motifj)

            motifs = motif

            if hm.score_motifs(motifs) < score_best_motifs:
                best_motifs = motifs
                score_best_motifs = hm.score_motifs(motifs)

        return best_motifs

    def create_profile(self, motifs):
        k = len(motifs[0])
        profile = {
            'A': [0]*k,
            'C': [0]*k,
            'G': [0]*k,
            'T': [0]*k,
        }

        for nuc in profile:
            for i in xrange(k):
                profile[nuc][i] = sum([1 if motif[i] == nuc else 0 for motif in motifs])*1.0/len(motifs)

        return profile

    def score_motifs(self, motifs):
        transposed = [''.join([motif[i] for motif in motifs]) for i in xrange(len(motifs[0]))]

        return sum([len(i) - Counter(i).most_common(1)[0][1] for i in transposed])


    def profile_most_probable_kmer(self, dna, k, pA, pC, pG, pT):
        profiles = {
            'A': pA,
            'C': pC,
            'G': pG,
            'T': pT
        }
        highest_prob = -1 
        highest_prob_kmer = ''

        for i in xrange(len(dna) - k + 1):
            current_kmer = dna[i:i+k]
            current_kmer_prob = reduce(lambda a, b: a*b, [profiles[current_kmer[nuc_i]][nuc_i] for nuc_i in xrange(len(current_kmer))])

            if current_kmer_prob > highest_prob:
                highest_prob = current_kmer_prob
                highest_prob_kmer = current_kmer

        return highest_prob_kmer

    def all_median_strings(self, dna, k):
        distance = 9999999999
        median = []

        for pattern in itertools.product('ACTG', repeat=k):
            pattern = ''.join(pattern)
            d_pattern_dna = sum([min([self.hamming_distance(pattern, dna_i[j:j+k]) for j in range(len(dna_i) - k + 1)]) for dna_i in dna])

            if distance == d_pattern_dna:
                median.append(pattern)
            elif distance > d_pattern_dna:
                distance = d_pattern_dna
                median = [pattern]

        return median

    def median_string(self, dna, k):
        distance = 9999999999
        median = ''

        for pattern in itertools.product('ACTG', repeat=k):
            pattern = ''.join(pattern)
            d_pattern_dna = sum([min([self.hamming_distance(pattern, dna_i[j:j+k]) for j in range(len(dna_i) - k + 1)]) for dna_i in dna])

            if distance > d_pattern_dna:
                distance = d_pattern_dna
                median = pattern

        return median

    def motif_enumeration(self, dna_set, k, d):
        patterns = set()
        raw_patterns = set()
        for dna in dna_set:
            for raw_pattern in [dna[i:i+k] for i in range(len(dna) - k + 1)]:
                raw_patterns.add(raw_pattern)

        for pattern in raw_patterns:
            for pattern_prime in self.neighbors(pattern, d):
                appears = True
                for dna in dna_set:
                    if self.approx_pattern_count(pattern_prime, dna, d) == 0:
                        appears = False
                        break

                if appears:
                    patterns.add(pattern_prime)

        return patterns

    # From: http://stackoverflow.com/a/19941659
    def neighbors(self, word, hamming_distance, charset='ATCG'):
        for indices in itertools.combinations(range(len(word)), hamming_distance):
            for replacements in itertools.product(charset, repeat=hamming_distance):
                mutation = list(word)
                for index, replacement in zip(indices, replacements):
                    mutation[index] = replacement
                yield "".join(mutation)

    # This is the brute force approach
    def freq_words_w_mismatches_and_reverse_comp(self, text, k, d):
        freq = [0] * (4**k)

        # Calculate freq
        for i in xrange(4**k):
            pattern = self.number_to_pattern(i, k)
            freq[i] = self.approx_pattern_count(pattern, text, d) + self.approx_pattern_count(self.reverse_complement(pattern), text, d)

        # Return only max words
        maxFreq = max(freq)
        maxWords = []
        for i in xrange(4**k):
            if freq[i] == maxFreq:
                maxWords.append(self.number_to_pattern(i, k))

        return maxWords

    # This is the brute force approach
    def freq_words_w_mismatches(self, text, k , d):
        freq = [0] * (4**k)

        # Calculate freq
        for i in xrange(4**k):
            freq[i] = self.approx_pattern_count(self.number_to_pattern(i, k), text, d)

        # Return only max words
        maxFreq = max(freq)
        maxWords = []
        for i in xrange(4**k):
            if freq[i] == maxFreq:
                maxWords.append(self.number_to_pattern(i, k))

        return maxWords

    def approx_pattern_count(self, pattern, text, d):
        return len(self.approx_pattern_matching(pattern, text, d))

    def approx_pattern_matching(self, pattern, text, d):
        len_pattern = len(pattern)
        pos = []
        for i in xrange(len(text) - len_pattern + 1):
            if self.hamming_distance(text[i:i+len_pattern], pattern) <= d:
                pos.append(i)

        return pos

    def hamming_distance(self, p, q):
        hamming_distance = 0
        for i in xrange(len(p)):
            if p[i] != q[i]:
                hamming_distance += 1

        return hamming_distance

    def min_skew(self, genome):
        min_skew_value = 0
        min_skew_i = []
        skew_array = [0]

        for i in xrange(len(genome)):
            base = genome[i]
            if base == 'C':
                skew_array.append(skew_array[-1] - 1)
            elif base == 'G':
                skew_array.append(skew_array[-1] + 1)
            else:
                skew_array.append(skew_array[-1])

            if skew_array[-1] < min_skew_value:
                min_skew_value = skew_array[-1]
                min_skew_i = [i + 1]
            elif skew_array[-1] == min_skew_value:
                min_skew_i.append(i + 1)

        return min_skew_i 

    def skew_genome(self, genome):
        skew_array = [0]

        for base in genome:
            if base == 'C':
                skew_array.append(skew_array[-1] - 1)
            elif base == 'G':
                skew_array.append(skew_array[-1] + 1)
            else:
                skew_array.append(skew_array[-1])

        return skew_array

    def pattern_to_number2(self, pattern):
        if not pattern:
            return 0

        symbol = pattern[-1]
        prefix = pattern[:-1]
        return 4 * self.pattern_to_number2(prefix) + self.lex_order[symbol]

    def computing_frequencies(self, text, k):
        freq = [0] * (4**k)
        for i in xrange(len(text) - k + 1):
            pattern = text[i:i+k]
            j = self.pattern_to_number(pattern)
            #print '%s -> %s'%(pattern, j)
            freq[j] += 1

        return freq

    def number_to_pattern(self, index, k):
        # 11, 2 = GT
        def recurse(num):
            if num < 4:
                return str(num)
            else:
                return str(recurse(num/4)) + str(num%4)

        base4_in_dna = [self.reverse_lex_order[i] for i in recurse(index)]
        return ''.join(['A']*(k-len(base4_in_dna)) + base4_in_dna)
        

    def pattern_to_number(self, pattern):
        # This is essentially base 4
        # GT = 11
        number = 0
        for i in xrange(len(pattern)):
            number += 4**(len(pattern) - 1 - i) * self.lex_order[pattern[i]]

        return number

    def clump_finder2(self, genome, k, L, t):
        # Generate freq array 
        freq_array = {}

        for i in range(len(genome) - k + 1):
            if genome[i:i+k] not in freq_array:
                freq_array[genome[i:i+k]] = [i]
            else:
                freq_array[genome[i:i+k]].append(i)

        # Calculate Lt-clumps
        tWords = []
        for word in freq_array:
            if len(freq_array[word]) < t:
                continue
            else:
                for i in range(len(freq_array[word]) - t + 1):
                    if freq_array[word][i+t-1] + k - freq_array[word][i] <= L:
                        tWords.append(word)
                        break

        return tWords
        

    def clump_finder(self, genome, k, L, t):
        def frequent_words4(text, k, t):
            # Do everything in one pass
            wordCounts = {}
            tWords = {};
        
            for i in range(len(text) - k + 1):
                currWord = text[i:i+k]
                
                # Update freq count of word in currWords dict
                if currWord in wordCounts:
                    wordCounts[currWord] += 1
                else:
                    wordCounts[currWord] = 1
        
                # See whether this is a new tWord entry
                if wordCounts[currWord] >= t:
                    tWords[currWord] = True;
        
            return tWords.keys()

        kLt_patterns = []
        for i in xrange(len(genome) - L):
            window_genome = genome[i:i+L]
            new_patterns = frequent_words4(window_genome, k, t)
            kLt_patterns += new_patterns

        return set(kLt_patterns)

    def pattern_matching(self, pattern, genome):
        patternLength = len(pattern)
        positions = []
    
        for i in range(len(genome)):
            if genome[i:i+patternLength] == pattern:
                positions.append(i)
    
        return ' '.join([str(i) for i in positions])

    def reverse_complement(self, pattern):
        return ''.join([self.substitution[i] for i in pattern][::-1])

    def frequent_words(self, text, k):
        # Generate all words of len k from text
        words = {}
        for i in range(len(text) - k + 1):
            words[text[i:i+k]] = 0
    
        # Count freq of words
        # and add entry into wordsByFreq as well
        wordsByFreq = {}
        for word in words.keys():
            words[word] = self.pattern_count(text, word)
    
            if words[word] not in wordsByFreq:
                wordsByFreq[words[word]] = [word]
            else:
                wordsByFreq[words[word]].append(word)
    
        # Return maximum
        return ' '.join(wordsByFreq[max(wordsByFreq.keys())])
    
    def frequent_words2(self, text, k):
        # Do everything in one pass
        wordCounts = {}
    
        for i in range(len(text) - k + 1):
            currWord = text[i:i+k]
            if currWord in wordCounts:
                wordCounts[currWord] += 1
            else:
                wordCounts[currWord] = 1
    
        maxFreq = max(wordCounts.values())
        return ' '.join([i for i in wordCounts.keys() if wordCounts[i] == maxFreq])
    
    def frequent_words3(self, text, k):
        # Do everything in one pass
        wordCounts = {}
        currMax = 0;
        maxWords = [];
    
        for i in range(len(text) - k + 1):
            currWord = text[i:i+k]
    
            # Update freq count of word in currWords dict
            if currWord in wordCounts:
                wordCounts[currWord] += 1
            else:
                wordCounts[currWord] = 1
    
            # See whether this is a new maxWord entry
            if wordCounts[currWord] == currMax:
                maxWords.append(currWord)
            elif wordCounts[currWord] > currMax:
                currMax = wordCounts[currWord]
                maxWords = [currWord]
    
        return ' '.join(maxWords)

    def pattern_count(self, text, pattern):
        patternLength = len(pattern)
        count = 0
    
        for i in range(len(text)):
            if text[i:i+patternLength] == pattern:
                count += 1
    
        return count

if __name__ == "__main__":
    hm = HiddenMessages()

    parser = argparse.ArgumentParser(description='Coursera Hidden Messages - Jhalley')
    parser.add_argument('--pattern_count', help='Count the number of occurances of pattern in text')
    parser.add_argument('--freq_words', help='Find most frequent k-mers in text')
    parser.add_argument('--freq_words2', help='Find most frequent k-mers in text - faster')
    parser.add_argument('--freq_words3', help='Find most frequent k-mers in text - fastest!')
    parser.add_argument('--reverse_complement', help='Find the reverse complement of a DNA string')
    parser.add_argument('--pattern_matching', help='Find all occurrences of a pattern in a string')
    parser.add_argument('--clump_finder', help='Find patterns forming clumps in a string')
    parser.add_argument('--clump_finder2', help='Find patterns forming clumps in a string using freq array')
    parser.add_argument('--computing_frequencies', help='Generate frequency array')
    parser.add_argument('--pattern_to_number2', help='Pattern to number algorithm')
    parser.add_argument('--number_to_pattern', help='My number to pattern algorithm')
    parser.add_argument('--skew_genome', help='Get skew diagram of genome')
    parser.add_argument('--min_skew', help='Find a position in a genome where the skew diagram attains a minimum')
    parser.add_argument('--hamming_distance', help='Calculate hamming distance between two strings of equal length')
    parser.add_argument('--approx_pattern_matching', help='Find all approximate occurrences of a pattern in a string')
    parser.add_argument('--approx_pattern_count', help='Get count of all approximate occurrences of a pattern in a string')
    parser.add_argument('--freq_words_w_mismatches', help='Find the most frequent k-mers with mismatches in a string')
    parser.add_argument('--freq_words_w_mismatches_and_reverse_comp', help='Find the most frequent k-mers (with mismatches and reverse complements) in a string')
    parser.add_argument('--motif_enumeration', help='A brute force attempt at motif finding')
    parser.add_argument('--median_string', help='Find a median string')
    parser.add_argument('--profile_most_probable_kmer', help='Find a Profile-most probable k-mer in a string')
    parser.add_argument('--greedy_motif_search', help='Greedy motif search')
    parser.add_argument('--greedy_motif_search_w_pseudocounts', help='Greedy motif search with pseudocounts')
    parser.add_argument('--randomized_motif_search', help='Randomized motif search')
    parser.add_argument('--gibbs_sampler', help='Gibbs sampler')
    args = parser.parse_args()

     
    # Parse from file
    if args.pattern_count:
        with open(args.pattern_count) as f:
            lines = f.readlines()
        print hm.pattern_count(lines[0].strip(), lines[1].strip())
    elif args.freq_words:
        with open(args.freq_words) as f:
            lines = f.readlines()
        print hm.frequent_words(lines[0].strip(), int(lines[1].strip()))
    elif args.freq_words2:
        with open(args.freq_words2) as f:
            lines = f.readlines()
        print hm.frequent_words2(lines[0].strip(), int(lines[1].strip()))
    elif args.freq_words3:
        with open(args.freq_words3) as f:
            lines = f.readlines()
        print hm.frequent_words3(lines[0].strip(), int(lines[1].strip()))
    elif args.reverse_complement:
        with open(args.reverse_complement) as f:
            lines = f.readlines()
        print hm.reverse_complement(lines[0].strip())
    elif args.pattern_matching:
        with open(args.pattern_matching) as f:
            lines = f.readlines()
        print hm.pattern_matching(lines[0].strip(), lines[1].strip())
    elif args.clump_finder:
        with open(args.clump_finder) as f:
            lines = f.readlines()
        klt = [int(i) for i in lines[1].strip().split()]
        print ' '.join(hm.clump_finder(lines[0].strip(), klt[0], klt[1], klt[2]))
    elif args.clump_finder2:
        with open(args.clump_finder2) as f:
            lines = f.readlines()
        klt = [int(i) for i in lines[1].strip().split()]
        print ' '.join(hm.clump_finder2(lines[0].strip(), klt[0], klt[1], klt[2]))
    elif args.computing_frequencies:
        with open(args.computing_frequencies) as f:
            lines = f.readlines()
        print ' '.join([str(i) for i in hm.computing_frequencies(lines[0].strip(), int(lines[1].strip()))])
    elif args.pattern_to_number2:
        with open(args.pattern_to_number2) as f:
            lines = f.readlines()
        print hm.pattern_to_number2(lines[0].strip())
    elif args.number_to_pattern:
        with open(args.number_to_pattern) as f:
            lines = f.readlines()
        print hm.number_to_pattern(int(lines[0].strip()), int(lines[1].strip()))
    elif args.skew_genome:
        with open(args.skew_genome) as f:
            lines = f.readlines()
        print hm.skew_genome(lines[0].strip())
    elif args.min_skew:
        with open(args.min_skew) as f:
            lines = f.readlines()
        print ' '.join([str(i) for i in hm.min_skew(lines[0].strip())])
    elif args.hamming_distance:
        with open(args.hamming_distance) as f:
            lines = f.readlines()
        print hm.hamming_distance(lines[0].strip(), lines[1].strip())
    elif args.approx_pattern_matching:
        with open(args.approx_pattern_matching) as f:
            lines = f.readlines()
        print ' '.join([str(i) for i in hm.approx_pattern_matching(lines[0].strip(), lines[1].strip(), int(lines[2].strip()))])
    elif args.approx_pattern_count:
        with open(args.approx_pattern_count) as f:
            lines = f.readlines()
        print hm.approx_pattern_count(lines[0].strip(), lines[1].strip(), int(lines[2].strip()))
    elif args.freq_words_w_mismatches:
        with open(args.freq_words_w_mismatches) as f:
            lines = f.readlines()
        print ' '.join(hm.freq_words_w_mismatches(lines[0].strip(), int(lines[1].strip().split()[0]), int(lines[1].strip().split()[1])))
    elif args.freq_words_w_mismatches_and_reverse_comp:
        with open(args.freq_words_w_mismatches_and_reverse_comp) as f:
            lines = f.readlines()
        print ' '.join(hm.freq_words_w_mismatches_and_reverse_comp(lines[0].strip(), int(lines[1].strip().split()[0]), int(lines[1].strip().split()[1])))
    elif args.motif_enumeration:
        with open(args.motif_enumeration) as f:
            lines = f.readlines()
        print ' '.join(sorted(hm.motif_enumeration([dna.strip() for dna in lines[1:]], int(lines[0].split()[0]), int(lines[0].split()[1]))))
    elif args.median_string:
        with open(args.median_string) as f:
            lines = f.readlines()
        print hm.median_string([dna.strip() for dna in lines[1:]], int(lines[0].strip()))
    elif args.profile_most_probable_kmer:
        with open(args.profile_most_probable_kmer) as f:
            lines = f.readlines()
        print hm.profile_most_probable_kmer(lines[0].strip(), int(lines[1].strip()), [float(i) for i in lines[2].strip().split()], [float(i) for i in lines[3].strip().split()], [float(i) for i in lines[4].strip().split()], [float(i) for i in lines[5].strip().split()])
    elif args.greedy_motif_search:
        with open(args.greedy_motif_search) as f:
            lines = f.readlines()
        print '\n'.join(hm.greedy_motif_search([dna.strip() for dna in lines[1:]], int(lines[0].strip().split()[0]), int(lines[0].strip().split()[1])))
    elif args.greedy_motif_search_w_pseudocounts:
        with open(args.greedy_motif_search_w_pseudocounts) as f:
            lines = f.readlines()
        print '\n'.join(hm.greedy_motif_search_w_pseudocounts([dna.strip() for dna in lines[1:]], int(lines[0].strip().split()[0]), int(lines[0].strip().split()[1])))
    elif args.randomized_motif_search:
        with open(args.randomized_motif_search) as f:
            lines = f.readlines()
        print '\n'.join(hm.randomized_motif_search(int(lines[0].strip().split()[0]), int(lines[0].strip().split()[1]), [dna.strip() for dna in lines[1:]]))
    elif args.gibbs_sampler:
        with open(args.gibbs_sampler) as f:
            lines = f.readlines()
        print '\n'.join(hm.gibbs_sampler(int(lines[0].strip().split()[0]), int(lines[0].strip().split()[1]), int(lines[0].strip().split()[2]), [dna.strip() for dna in lines[1:]]))

    # Test calls
    # print hm.score_motifs(['TCTCGGGG', 'CCAAGGTG', 'TACAGGCG', 'TTCAGGTG', 'TCCACGTG'])
    print '\n'.join(hm.gibbs_sampler(8, 5, 100, ['CGCCCCTCTCGGGGGTGTTCAGTAACCGGCCA', 'GGGCGAGGTATGTGTAAGTGCCAAGGTGCCAG', 'TAGTACCGAGACCGAAAGAAGTATACAGGCGT', 'TAGATCAAGTTTCAGGTGCACGTCGGTGAACC', 'AATCCACCAGCTCCACGTGCAATGTTGGCCTA']))
    #print Counter([hm.random_biased_die([0.1, 0.2, 0.3]) for i in xrange(10000)])
    #print hm.randomized_motif_search(8, 5, ['CGCCCCTCTCGGGGGTGTTCAGTAAACGGCCA', 'GGGCGAGGTATGTGTAAGTGCCAAGGTGCCAG', 'TAGTACCGAGACCGAAAGAAGTATACAGGCGT', 'TAGATCAAGTTTCAGGTGCACGTCGGTGAACC', 'AATCCACCAGCTCCACGTGCAATGTTGGCCTA'])
    # print hm.all_median_strings(['CTCGATGAGTAGGAAAGTAGTTTCACTGGGCGAACCACCCCGGCGCTAATCCTAGTGCCC', 'GCAATCCTACCCGAGGCCACATATCAGTAGGAACTAGAACCACCACGGGTGGCTAGTTTC', 'GGTGTTGAACCACGGGGTTAGTTTCATCTATTGTAGGAATCGGCTTCAAATCCTACACAG'], 7)
{'fraction_non_alphanumeric': 0.12247947641206068, 'fraction_numerical': 0.010781977074111906, 'mean_word_length': 2.1841835812804913, 'pattern_counts': {'":': 2, '<': 28, '<?xml version=': 0, '>': 17, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 4, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '6768302', 'n_tokens_mistral': 18073, 'n_tokens_neox': 17074, 'n_words': 3825}