# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3.9.0 64-bit
#     name: python3
# ---

# ## General Notes on Motif Selection for Transcription Factor of Circadian Rythm
#
# - Circadian clock must have some basis to regulate body functions → this can be due to a clock gene which regulates gene expression and consequently protein synthesis for certain bodily effects → we see similar correlations where some heart attacks happen more often in morning and asthma attacks more common at night
# - 1970 <NAME> and <NAME> identified mutant flies with abnormal circadian patterns and traced fly mutations to a single regulatory gene or transcription factor → 2 decades after and beyond, more circadian genes have been discovered in mammals
# - these specific genes have shown high degree of evolutionary conservation across species
# - every plant cell keeps track of day and night independently → regulatory proteins listed initially are controlled by external factors like nutrient availability or sunlight → and allow for individual cells to specifically adjust certain gene expression
#     - TOC1 promotes expression of LHY and CCA1, LHY and CCA1 repress expression of TOC1 → this is a **negative feedback loop**
#     - sunlight triggers transcription of LHY and CCA1, thus represses TOC1 transcription
#     - based on this, we know that transcription of TOC1 peaks at midnight and slowly starts to promote the expression of CCA1 and LHY which can only be promote once sunlight increases
# - the reason these 3 genes can control transcription and expression of other genes because they are transcription factors known as regulatory proteins that turn other genes on and off
#     - regulates gene by binding to specific short DNA interval called a **regular motif also known as a transcription factor binding site** **in the gene's upstream region** → 600-1000 nucleotide-long region preceeding start of the gene
#     - CCA1 binds to AAAAAATCT in upstream region of many genes regulated by CCA1
#     - these motifs **are not consevered** and subject to mutations naturally → positions also vary primarly
#         - CCA1 may bind to AAGAACTCT instead of the traditional binding site
#
# **How can we locate these regulatory motifs/binding sites without knowing what they look at in advance?** 
#
# - you can easily validate this experimentally, but that takes time and resources and automation and pattern finding could be a much easier approach to finding this → how can we find **AAAAAATCT for CCA1 with just bioinformatics**
#
# **THE MOTIF FINDING PROBLEM:** 
#
# *Experimentation:* <NAME> uses DnA arrays to see what genes are activated based on the time of day
#
# - did this by extracting around 500 genes that exhibited changes based on circadian behaviour and looked for the most frequent patterns in the gene's upstream regions (***upstream = end of coding/3'→5' direction)**
# - found that AAAATATCT is frequent motif appearing 46 times
#
# For each of the 500 genes, he concatenated the upstream regions which actually is not good because it reduces the potential for accurate pattern finding between each window of the region + it reduces the chance that there could be mutations → easier to find explicit patterns than nuances
#
# **Concatenated Length →** 500 regions each 1000 nucleotides long, 9-mer motif
#
# - given there are 4 possible nucleotides and 9 length, the total permutation of possible 9-mers is 4^9 = 262144
# - likelihood of our most favourable 9-mer is 1/262144 →
# - per region, we can get 1000-(9-1) possible 9-mers which is 992, and this applies for the all 500
# - therefore, total probability is 500 * 992 * 1/262144 = **1.8921**
#
# - to validate his finding, Kay mutated the sequence in the upstream and saw that the gene lost its circadian behaviour
# - the issue is that most motifs are elsuze → fly infected with bacterium develops immunity gene expression, where some genes have elevated expression levels and are likely to be immunity genes → they bind to the transcription factor NF-kB but these sites are not conserved or the same throughout → the following shows ten NF-kB binding sites from Drosophila melanogaster genome
#
#     <img src='https://s3.us-west-2.amazonaws.com/secure.notion-static.com/917b2b8a-2675-4f6c-a7bd-f3f2d485545e/Screen_Shot_2021-06-15_at_12.48.52.png?X-Amz-Algorithm=AWS4-HMAC-SHA256&X-Amz-Credential=AKIAT73L2G45O3KS52Y5%2F20210615%2Fus-west-2%2Fs3%2Faws4_request&X-Amz-Date=20210615T165800Z&X-Amz-Expires=86400&X-Amz-Signature=db03516ca0af99358cdd01564d0e81585f06d23908194a9af526f76741198422&X-Amz-SignedHeaders=host&response-content-disposition=filename%20%3D%22Screen_Shot_2021-06-15_at_12.48.52.png%22' />
#
#     - there is a lot of variation in these binding sites, and because there are so many possible outcomes, it becomes difficult to identify one because we don't have an initial/template column
#         - ie. if we chose one of these as a template and constrainted the conserved template to the nucleotides that are common for each and every factor, the mutation margin error will prevent us from finding all possible NF-kB factors
#     - **WE CANNOT USE FREQUENT WORDS** because variation of the regulatory motif is not a distinctive trait as it can appear sparsely over just a handful of regions scattered across the genome

# +
##Brute Force Implementation: 

import numpy as np

def hamming_distance(p, q):
    value_dict = {'A':1, 'T':2, 'C':3 ,'G':4}
    new_p = np.array([value_dict[k] for k in p])
    new_q = np.array([value_dict[k] for k in q])
    
    diff = np.subtract(new_p, new_q)
    return len(diff) - list(diff).count(0)

def mutate_sequence(sequence, k):
    mutated_sequence = []
    nucleotides = ['A', 'C', 'T', 'G']
    for idx, nuc in enumerate(sequence):
        all_possibles = []
        all_possibles += [sequence[:idx] + i + sequence[idx+1:] for i in nucleotides]
        mutated_sequence+= all_possibles
    return mutated_sequence

def generate_neighbours(sequence, k):
    all_neighbours = []
    for i in range( 0, len(sequence)-k+1):
        base_neighbours = sequence[i:i+k]
        all_neighbours.append(base_neighbours)
    return all_neighbours

def MotifEnumeration(sequences, k, d):
    all_patterns = {item: set() for item in sequences}
    
    for idx, sequence in enumerate(sequences):
        base_neighbours = generate_neighbours(sequence, k)
        mutated_neighbours = []
        for neighbour in base_neighbours:
            mutated_neighbours += mutate_sequence(neighbour, k)
        all_patterns[sequence] = set(mutated_neighbours)
        
    
    return all_patterns

# +
with open('dataset_156_8.txt', 'r') as sequence:
    txtinput = sequence.read().splitlines()
    
params = [int(i) for i in txtinput[0].split(' ')]
sqinput = txtinput[1:]


all_neighbours = MotifEnumeration(sqinput, params[0], params[1])
patterns = list(all_neighbours.values())

base_patterns = patterns[0]
for i in patterns[1:]:
    base_patterns.intersection_update(i)
print(*base_patterns)
# -

# ## The inefficiences of the algorithm above are clearly unscalable. We should solve another problem
# - implanted motif problem is a useful abstraction of biological problem of motif finding, has limitations 
# - during experimentation, the fair assumption was made that NF-kB is not the only transcription binding region that directly regulates gene expression
#     - some genes can be influenced by more than just the circadian rhythm
# - scoring motifs requires a list of t DNA strings of Dna, where you select a k-mer from each string to form a collection Motifs -> (t * k) motif matrix
# - we score motifs by counting the total number of nucleotides that do not match the most frequent nucleotide in said column, after which you can sum up the scores
# - we count mofifs by looking at the frequency of each nucleotide in the column 
#     - you can use the count and divide it by the number of rows to gain a percentage point decimal of that given nucleotide and its makeup of the entire column for all potentian transcription sites
#     
# - based on certain observations, most biological motifs have 2 nucleotides that have a conservative position giving them the ability to bind to a transcription factor
#     - 16 nucleotide-long CSRE transcription factor binding site consists of five positions and 11 weak positions
#     ![Screen Shot 2021-06-18 at 15.46.43.png](attachment:f30c9248-96a0-4a6e-a6b5-2a6eb1c1ea06.png)
#     - when comparing two columns, we want to have mismatched nucleotides that are better alternatives/more viable for our conservative strand
# - therefore, (6C, 4T) is more conserved than (6C, 2A, 2T) and should recieve a lower score -> another way to think of this is variation of our probability distribution
#
# - model the scroe as probability distribution -> ![Screen Shot 2021-06-18 at 15.50.31.png](attachment:c45e90a4-3606-4993-a6b5-0cf03bf4744a.png)
#     - to calculate variation or the *entropy* of the probability distribution, use ![Screen Shot 2021-06-18 at 15.54.22.png](attachment:f5d7f944-0a1a-45a0-9c95-3496efe07f63.png)
#         - the formula is the negative sum of the products of a given probability distribution profile and its log-base-2 -> this is ultimately showing just how uncertain the chances are
#     - the lower the entropy, the less uncertainty and variation there is in the distribution, or the more conservative the strand is
#

# +
# Entropy Function
import math


def Count(motif_matrix, laplace=True): 
    motif_matrix = list([i for i in motif_matrix])
    N = len(motif_matrix[0]) * 2
    count_matrix = [[] for _ in range(len(motif_matrix))]
    if laplace == True:
        for idx, column in enumerate(motif_matrix):
            count_matrix[idx] = np.asarray([np.asarray((list(column).count(i)+1)/N) for i in ['A', 'C', 'G', 'T']])
    elif laplace == False:
        N = len(motif_matrix[0])
        for idx, column in enumerate(motif_matrix):
            count_matrix[idx] = np.asarray([np.asarray((list(column).count(i))/N) for i in ['A', 'C', 'G', 'T']])
    return np.asarray(count_matrix)

def H_e(counts):
    base = np.nan_to_num(counts*np.log2((counts)))
    sums = [-1 * sum(base[i]) for i in range(len(base))]
    return sums


# +
lst = ['TCGGGGGTTTTT', 'CCGGTGACTTAC', 'ACGGGGATTTTC', 'TTGGGGACTTTT', 'AAGGGGACTTCC', 'TTGGGGACTTCC', 'TCGGGGATTCAT', 'TCGGGGATTCCT', 'TAGGGGAACTAC', 'TCGGGTATAACC']

formatted = [np.array([i for i in x]) for x in lst]
formatted = np.array([i for i in formatted]).transpose()

counts = Count(formatted, laplace=False)
print(counts.transpose())
entropy = H_e(counts)
# print(entropy)
print(sum(entropy))


# -

# ### The Motif Finding Probelm
# Given a collection of strings, find a set of kmers that minimzes the score of the resulting motif -> Inp = collections of strings and int-k  |  Output = collection of motifs of k-mers, minimizing Score(Motifs)
#
# Because this is a brute-force approach, we should redefine the problem into something simpler (use hamming distance function to optimize our approach by degenerate sequence with consensus sequence). 
#
# d(Pattern, Motif(i-n)) = sum(HammingDistance(Pattern, Motif(i))
#
# <strong>Instead of generating a random set of motifs that assume the lowest score, we should try to calculate the hamming score of the pattern/selected k-mer that minimizes that hamming score across the entire sequence.</strong>
#

# +

def AllKmers(k):
    import itertools
    possibles = 4**k
    init = 'ACTG'
    all_sequences = set(map("".join, itertools.product(init, repeat=k)))
    return all_sequences
    

def MedianString(k, dna_list):
    min_pattern = {}
    
    all_sequences = AllKmers(k)
    for i in all_sequences:
        pattern_score = []
        for x in dna_list:
            all_neighbours = generate_neighbours(x, k)
            median_score = min([hamming_distance(i, z) for z in all_neighbours])
            pattern_score.append(median_score)
        min_pattern[i] = sum(pattern_score)
    
    return min_pattern
            


# +
with open('dataset_158_9.txt', 'r') as seqlist:
    seqinp = seqlist.read().splitlines()
    k_len = int(seqinp[0])
    seqinp = seqinp[1:]

pattern_min = MedianString(k_len, seqinp)
print(pattern_min)

# +
min_value = min(pattern_min.values())
vals = list(pattern_min.values())
keys = list(pattern_min.keys())

print(keys[vals.index(min_value)])


# -

# ## This stuff is too slow, we should move onto Greedy Motif Search
# Most Brute-force algorithms iteratively generate possibilities and check whether they satisfy a condition, but this is inefficient design of algorithms and almost a shortcut. Therefore, greedy algorithms are favourable as they select the most fit solution by checking the criteria at each iteration. 
#
# This criteria helps the algorithm find a direction/make a choice at each iteration to optimize its solution to the problem. This is why greedy algorithms often make a locally optimum solution. At each iteration, it gets easier and easier to make decisions (if not any changes at all). 
#
# Most applications can only be done with optimal substructure problems where there will always be an optimal solution found in a simplified problem along with the local minimum typically being cloes if not equal to the global optimum. 
#

# +
def ProfileProbable(text, k, profile):
    all_neighbours = list(set(generate_neighbours(text, k)))
    profile_map = dict(zip(['A', 'C', 'G', 'T'], profile))
    prod_vals = [1 for _ in range(len(all_neighbours))]
    indices = [i for i in range(0, k)]
    motif_profile = {motif: [] for motif in all_neighbours}
    
    
    for motif in motif_profile.keys():
        template = []
        prod = 1
        for idx, nuc in enumerate([x for x in motif]):
            template.append(profile_map[str(nuc)][int(idx)])
        for i in template:
            prod *= float(i)
        motif_profile[motif] = prod
    
    return motif_profile, profile_map, list(motif_profile.keys())[list(motif_profile.values()).index(max(list(motif_profile.values())))]
    

with open('dataset_159_3.txt', 'r') as sqinput:
    sqinput = sqinput.read().splitlines()
    text = sqinput[0]
    k = int(sqinput[1])
    profile_matrix = [sqinput[i].split(' ') for i in range(2, len(sqinput))]
    
# profile_matrix = ['0.2 0.2 0.3 0.2 0.3'.split(' '), 
#                   '0.4 0.3 0.1 0.5 0.1'.split(' '), 
#                   '0.3 0.3 0.5 0.2 0.4'.split(' '), 
#                   '0.1 0.2 0.1 0.1 0.2'.split(' '), ]

# for i, row in enumerate(profile_matrix):
#     for x, item in enumerate(row):
#         profile_matrix[i][x] = float(item)

entropy, profile_map, probable = ProfileProbable(text, k, profile_matrix)


# +
def NewProfileProbable(text, k, profile):
    all_neighbours = generate_neighbours(text, k)
    profile_map = dict(zip(['A', 'C', 'G', 'T'], profile))
    all_scores = {motif: 0 for motif in all_neighbours}
    
    for motif in all_neighbours:
#         print(motif)
        score_motif = []
        
        for idx, nuc in enumerate(str(motif)):
            score_motif.append(profile_map[nuc][idx])
        
        prod = 1
        for i in score_motif:
            prod *= float(i)
        all_scores[motif] = prod
        
    return profile_map, list(all_scores.keys())[list(all_scores.values()).index(max(list(all_scores.values())))]


def GreedyMotifSearch(k, t, text):
#     profile_map = dict(zip(['A', 'C', 'G', 'T'], [[0, 0, 0], [0, 0, 0], [0, 0, 0], [0, 0, 0]]))
    gen_patterns = generate_neighbours(text[0], k)
    best_motifs = {key: [[]] for key in gen_patterns}
    all_scores = []
    count = 0
    for current_motif in list(best_motifs.keys()):  
        lst1 = [current_motif]
        for line in text[1:]:
            line = str(line)
            formatted = [np.array([i for i in x]) for x in lst1]
            formatted = list(np.array([i for i in formatted]).transpose())
            if lst == [current_motif]:
                profile_matrix = list(Count(formatted, laplace=True).transpose())
            else:
                profile_matrix = list(Count(formatted, laplace=False).transpose())e=False).transpose())
#             print(profile_matrix)
#             prin
est = NewProfileProbable(line, k, profile_matrix)
#             prilst1profile_map)
            best_motifs[current_motif][0].append(best)
            lst.append(best)
        entropy = H_e(profile_matrix)
        best
        # print(count)_motifs[current_motif].append(sum(entropy))
        all_scores.append(sum(entropy))
    return best_motifs, min(all_scores)
    
    

# + tags=[]
#Implement Laplace

text = ['GGCGTTCAGGCA', 'AAGAATCAGTCA', 'CAAGGAGTTCGC', 'CACGTCAATCAC', 'CAATAATATTCG']

with open('dataset_160_9 (65).txt', 'r') as sqinput:
    sqinput = sqinput.read().splitlines()
    params = [int(i) for i in sqinput[0].split(' ')]
   #  text = sqinput[1:]    
    
best_motif_dict, best_score = GreedyMotifSearch3] 5,, text)

# print(best_motif_dict)
best_motifs = [[i[0] for i in list(best_motif_dict.values()) if i[1] == best_score], best_score]
best_key = [key for key in best_motif_dict.keys() if list(best_motif_dict[key])[1] == best_score]
all_scores = [i[1] for i in list(best_motif_dict.values())]
# print(best_motif_dict)

print(*best_key)
print(*best_motif_dict[best_key[0]][0], sep='\n')
print(best_score)
# -

# ### We have to recognize that Greedy Motif Finding itself is flawed in this process primarily because the profile matrix differs when picking the most relevant motifs. 
#
# This means the algorithm can miss certain patterns because it requires much more data to formulate an accurate profile. When observing the second line in the sequence, anything that does not match the 1, 1, 0 profile matrix (assuming d=1) will be discarded but we all know this profile matrix is an incorrect representation of the data. And as you move along each line, the matrix changes and is not constant.
#
# Whenever we have a 0 probability, it makes the entire sequence's entropy 0. This skews and hides potentially important kmers. 
#
# We should try and utilize Laplace's Rule of Succession which adds 1 to each element of the total count of the nucleotides in a motif. For example, if your given probability was 10/10, Laplace's Rule of Succession makes you add 1 to the numerator/posterior and 2 to the denominator/sample_size. So now you go from 100% to 91.7%. But how does this play into our understanding of profile matrices and Greedy Search?
#
# To give some context, our goal is to look at a bunch of genes and their sequences and identify the motifs that most closely resemble *each other*. This is critical because we aren't looking for a most-frequent or anything specific like that. Instead we're just looking for how conserved motifs are collectively. The way we model this is by looking at all of our sequences, after which we generate a profile for the probability distribution of all possible motifs. What this does is that it gives us a good idea of the most frequent occuring nucleotide at a given position for our "ideal motif" as position and nucleotide matter in this case. Because there is a random chance of any mutation happening, it's hard to take the observed probability as the most accurate probability that accounts for a given error rate when new data is introduced. 
#
# In simulations, note that as more data is introduced, we get closer to the true error rate but that doesn't mean that it serves you well for a given case. 
# ![Screen Shot 2021-06-23 at 10.11.08.png](attachment:1120d861-bb9d-4d69-8717-06b33ec16dfa.png)
#
# To make this more effective, we only really need to initialize our profile matrix formation with an extra count. 
#
#

# +
# text = ['GGCGTTCAGGCA', 'AAGAATCAGTCA', 'CAAGGAGTTCGC', 'CACGTCAATCAC', 'CAATAATATTCG']

with open('dataset_160_9 (7).txt', 'r') as sqinput:
    sqinput = sqinput.read().splitlines()
    params = [int(i) for i in sqinput[0].split(' ')]
    text = sqinput[1:]    
    
best_motif_dict, best_score = GreedyMotifSearch(params[0], params[1], text)

# print(best_motif_dict)
best_motifs = [[i[0] for i in list(best_motif_dict.values()) if i[1] == best_score], best_score]
best_key = [key for key in best_motif_dict.keys() if list(best_motif_dict[key])[1] == best_score]
print(best_motif_dict)

print(*best_key)
print(best_score)

for x in [i for i in best_motifs[0]][0]:
    print(x)
print(best_motif_dict['GAAATCATACGC'])

# +
import math 

def D_ps(pattern, dna):
    k = len(pattern)
    distance = 0
    for line_text in dna:
        hamming_d = math.inf
        for motif in generate_neighbours(line_text, k):
            new_hamming_d = hamming_distance(pattern, motif)
            if hamming_d > new_hamming_d:
                hamming_d = new_hamming_d
        distance = distance + hamming_d
    return distance

D_ps('AAA', ['TTACCTTAAC', 'GATATCTGTC', 'ACGGCGTTCG', 'CCCTAAAGAG', 'CGTCAGAGGT'])

# +
with open('dataset_5164_1.txt', 'r') as sqinput:
    sqinput = sqinput.read().splitlines()
    pattern = sqinput[0]
    input_dna = sqinput[1].split(' ')

D_ps(pattern, input_dna)


# -

# ### Despite out best efforts, the code we previously had is flawed in very specific circumstances. There are 4 main errors that are the most likely cause of our failure: 
# 1. We aren't computing our score correctly for the overall sequence -> the correct score should be the entropy of the selected pattern and the sum of the respective pattern's entropys for a selected motif. 
# 2. The profile matrix is flawed. Perhaps the method of counting is flawed and the count -> probability is also wrong. This is most likely due to the fact that our denominator is wrong and initialization of the profile matrix is wrong. I noticed I'm only doing it once for the first selected pattern and not initializing it later. Additionally, it could be that the shape and profile matrix generated itself is wrong from the counts itself, so double checking that would be better. 
# 3. There is no flag to figure out issues when 2 or more pattern scores are the same. Ideally, we want to be able to print the first one that comes up so we need to write that down accordingly. 
#
# We are going to address this primarily by rewriting some functions and double checking if they have any inconsistencies with certain test cases. 

# +
##Because my greedy motif thing is not working that well, I'm going to try and rewrite the entire system here so wish me luck :(

def newneighbours(text, k):
    return [text[i:i+k] for i in range(len(text) - (k-1))]

def entropy_calc(counts):
    base = np.nan_to_num(counts*np.log2((counts)))
    sums = [-1 * sum(base[i]) for i in range(len(base))]
    return sums

def hamming_distance(p, q):
    value_dict = {'A':1, 'T':2, 'C':3 ,'G':4}
    new_p = np.array([value_dict[k] for k in p])
    new_q = np.array([value_dict[k] for k in q])
    
    diff = np.subtract(new_p, new_q)
    return len(diff) - list(diff).count(0)

def best_in_line(line, profile, k):
    min_ent = math.inf
    best_motif = ''
    all_possible_motifs = newneighbours(line, k)
    for motif in all_possible_motifs:
        score = 1
        for idx, nuc in enumerate(motif):
            print(nuc, idx, profile[nuc][idx])
        if min_ent > score:
            min_ent = score
            best_motif = motif

    return best_motif, min_ent

def greedy_cromwell(k, t, text):
    possible_motifs = newneighbours(text[0], k)
    best_selections = {p: [[], 1] for p in possible_motifs}

    for motif in possible_motifs:
        text_for_profile = [motif]
        base_entropy = 0
        for line in text[1:]:
            _, profile, _ = text_counts(text_for_profile, k, laplace=True)
            best_motif_line, entropy = best_in_line(line, profile, k)
            text_for_profile.append(best_motif_line)
            best_selections[motif][0].append(best_motif_line)
            base_entropy += entropy

        best_selections[motif][1] = base_entropy
    minima = min([i[1] for i in list(best_selections.values())])
    best_mot = [motif for motif in best_selections.keys() if best_selections[motif][1] == min([i[1] for i in list(best_selections.values())])]
    return best_selections, minima, best_mot


# +
DNA = ['GGCGTTCAGGCA', 'AAGAATCAGTCA', 'CAAGGAGTTCGC', 'CACGTCAATCAC', 'CAATAATATTCG']
k = 3
t = len(DNA)

best, mi, ma = greedy_cromwell(k, t, DNA)
print(best)
print(mi)
print(ma) #there are some primary errors in our profile matrix as seen below, we should try to analyze and trace-back our steps

# -

# ### To test my theories in errors, I decided to rewrite the first counting-profile method and then hard-code what the greedy algorithm would have done in the algorithm by following the example and logic previously discussed
#
# First step, rewrite the text_count function such that is accurately reflects the probability matrix we're hoping to generate. This was actually more useful than our previous one because we could easily change the calculation methodology for laplace-init, albeit it was slower due to the for-loops over the matrix representation.
#
# Through this process, I realized that a lot of my flaws came from generating the profile matrix and then picking the correct sequence of patterns as seen below. Through this loop, I found the correct process for picking the ideal motif and I implemented that later on in the newgreedyfind function. 

# +
def text_counts(text, k, laplace=True): #after rewriting the count and profile method, I noticed the main flaw in the profile_dict calculations
    if laplace == True:
        init = 1
    else:
        init = 0

    count_dict = {key: [init] * k for key in ['A', 'C', 'G', 'T']}

    for row in text:
        row = str(row)
        for idx, nuc in enumerate(row):
            count_dict[nuc][idx] += 1

    profile_dict = {key: [i/(k+(init*len(text))) for i in value] for key, value in count_dict.items()} #here, we need to divide by the len(text) or t + k as it matches the laplace-init
    di = list(count_dict.values())
    
    return count_dict, profile_dict, entropy_calc(np.asanyarray(di).transpose())

# DNA = ['GGCGTTCAGGCA', 'AAGAATCAGTCA', 'CAAGGAGTTCGC', 'CACGTCAATCAC', 'CAATAATATTCG']
# k = 3
# t = len(DNA)

DNA = ['TTACCTTAAC', 'GATGTCTGTC', 'ACGGCGTTAG', 'CCCTAACGAG', 'CGTCAGAGGT']
k = 4
t = len(DNA)



primary_list = ['ACCT']
count, profile, _ = text_counts(primary_list, k, laplace=True)
print(count)
print(profile)

scores_line_1 = {key: 0 for key in newneighbours(DNA[1], 4)}
max_prod = 0

for pattern in newneighbours(DNA[1], 4):
    prod = 1
    for idx, nuc in enumerate(pattern):
        prod *= profile[nuc][idx]
    scores_line_1[pattern] = prod
    
    if max_prod < prod:
        max_prod = prod

print(scores_line_1, max_prod)

first_motif = [key for key in scores_line_1.keys() if scores_line_1[key] == max_prod][0]
print(first_motif)
primary_list.append(first_motif)

count, profile, _ = text_counts(primary_list, k, laplace=True)
print(count)
print(profile)


scores_line_2 = {key: 0 for key in newneighbours(DNA[2], 4)}
max_prod = 0

for pattern in newneighbours(DNA[2], 4):
    prod = 1
    for idx, nuc in enumerate(pattern):
        prod *= profile[nuc][idx]
    scores_line_2[pattern] = prod
    
    if max_prod < prod:
        max_prod = prod

print(scores_line_2, max_prod)



# +

def newgreedyfind(k, t, text):
    init_motifs = newneighbours(text[0], k)
    primary_lists = {key: [[], []] for key in init_motifs}

    for pattern in init_motifs:
        primary_list = [pattern]
        for line in text[1:]:
            count, profile, entropy = text_counts(primary_list, k, laplace=True) 
            score_line = {key: 0 for key in newneighbours(line, k)}
            max_prod = 0
            for selected_pattern in newneighbours(line, k):
                prod = 1
                for idx, nuc in enumerate(selected_pattern):
                    prod *= profile[nuc][idx]
                score_line[selected_pattern] = prod 

                if max_prod < prod:
                    max_prod = prod
            
            first_motif = [key for key in score_line.keys() if score_line[key] == max_prod][0] #another key thing I missed was that we must generate the profile matrix from the best_motifs
            primary_lists[pattern][0].append(first_motif)
            primary_lists[pattern][1].append(max_prod)
            

    return primary_lists


with open('dataset_160_9 (8).txt', 'r') as sqinput:
    sqinput = sqinput.read().splitlines()
    params = sqinput[0].split(' ')

DNA = sqinput[1:]
k = int(params[0])
t = len(DNA)


pattern_best_selections = newgreedyfind(k, t, DNA)

best_sum_score = [0, 0]
best_prod_score = [0, 0]

all_sums = {key: 0 for key in pattern_best_selections.keys()}

for key, value in pattern_best_selections.items():
    pattern_best_selections[key][0] = pattern_best_selections[key][0][:t-1]
    pattern_best_selections[key][1] = pattern_best_selections[key][1][:t-1]

    all_sums[key] = sum(pattern_best_selections[key][1])


print(*[key for key in all_sums.keys() if all_sums[key] == max(all_sums.values())] + [pattern_best_selections[key][0] for key in all_sums.keys() if all_sums[key] == max(all_sums.values())][0])
    
# -



StarcoderdataJupyterScriptsDedupFiltered