asacas/CodonOpt
0
1from numpy.random import randint2from numpy.random import rand3import time4from scipy.optimize import minimize5from sklearn.model_selection import GridSearchCV6from sklearn.model_selection import RandomizedSearchCV7import pandas as pd8import numpy as np9import re10from collections import Counter11import math12import itertools13import pprint14import random15import copy 16import pygad17import statistics18 19def parse(file):20 file = re.sub(r'\( ', '', file)21 file = re.sub(r'\) ', '', file)22 file = re.sub(r'\n', '', file)23 file = re.sub(r'[(|)]', ' ', file)24 file_split = file.split(" ")25 26 rows = []27 row = []28 counter = 129 for i in file_split:30 #print(i, counter)31 if i == "":32 continue33 if counter == 5:34 #print(row)35 rows.append(row)36 row = []37 counter = 138 continue39 if counter < 5:40 row.append(i)41 counter += 142 43 profile = pd.DataFrame(columns = ('Codon', 'AA', 'Frequency', 'Usage'))44 profile = pd.DataFrame(rows).rename({0: "Codon", 1: "AA", 2: "Frequency", 3: "Usage"}, axis='columns')45 return(profile)46 47def read_and_parse(file_path):48 file = open(file_path, 'r').read()49 file = re.sub(r'\( ', '', file)50 file = re.sub(r'\) ', '', file)51 file = re.sub(r'\n', '', file)52 file = re.sub(r'[(|)]', ' ', file)53 file_split = file.split(" ")54 55 rows = []56 row = []57 counter = 158 for i in file_split:59 #print(i, counter)60 if i == "":61 continue62 if counter == 5:63 #print(row)64 rows.append(row)65 row = []66 counter = 167 continue68 if counter < 5:69 row.append(i)70 counter += 171 72 profile = pd.DataFrame(columns = ('Codon', 'AA', 'Frequency', 'Usage'))73 profile = pd.DataFrame(rows).rename({0: "Codon", 1: "AA", 2: "Frequency", 3: "Usage"}, axis='columns')74 return(profile)75 76def AA_look_up(AA_seq, profile):77 AA_seq = [x for x in AA_seq]78 AA_seq = pd.Series(AA_seq)79 80 counts = pd.DataFrame(AA_seq.value_counts())81 counts.index.name = 'AA'82 counts.reset_index(inplace = True)83 counts = counts.rename(columns = {0 : "Total"})84 look_up_df = pd.merge(counts, profile, left_on = "AA", right_on = "AA")85 look_up_df_type_dict = {86 'AA' : str,87 'Total': float,88 'Codon' : str,89 'Frequency': float,90 'Usage' : str91 }92 look_up_df = look_up_df.astype(look_up_df_type_dict)93 look_up_df['#'] = np.round(look_up_df['Total'] * look_up_df['Frequency'])94 look_up_df['Used'] = 095 96 look_up_arr = np.asarray(look_up_df)97 #print(np.unique(AA_seq))98 look_up_arr_list = []99 for AA in np.unique(AA_seq):100 sub_look_up_arr = look_up_arr[look_up_arr[:,0] == AA]101 #print(sub_look_up_arr)102 total = sub_look_up_arr[0,1]103 col_sum = sub_look_up_arr.sum(axis = 0)[5]104 #print(sub_look_up_arr, total, col_sum)105 #col_sum = col_sum - 5106 counter = total - col_sum107 if counter != 0:108 modifier = counter/abs(counter)109 #counter = 10110 #print(counter, modifier)111 #print(counter)112 #y = 113 freq = np.sort(sub_look_up_arr[:,3])114 freq[-1]115 if counter > 0:116 used_highest = sub_look_up_arr[sub_look_up_arr[:,5] == max(sub_look_up_arr[:,5]), 6] 117 sub_look_up_arr[sub_look_up_arr[:,5] == max(sub_look_up_arr[:,5]), 6] = used_highest + 1118 if counter < 0:119 used_highest = sub_look_up_arr[sub_look_up_arr[:,5] == min(sub_look_up_arr[:,5]), 6] 120 sub_look_up_arr[sub_look_up_arr[:,5] == min(sub_look_up_arr[:,5]), 6] = used_highest - 1 121 x = 0*len(sub_look_up_arr[:,6])122 #print(sub_look_up_arr[:,6] + 1, x)123 124 while counter != 0:125 for row in sub_look_up_arr:126 #print(row)127 if counter == 0:128 #print(sub_look_up_arr)129 break130 if row[6] != 0:131 continue132 133 134 sub_look_up_arr[sub_look_up_arr[:,2] == row[2], 5] += 1 * modifier135 sub_look_up_arr[sub_look_up_arr[:,2] == row[2], 6] += 1 * modifier136 137 #print(sub_look_up_arr, counter)138 139 if counter > 0:140 counter = counter - 1141 else:142 counter = counter + 1143 144 sub_look_up_arr[:,6] = 0145 sub_look_up_arr[:,6] = 100146 look_up_arr_list.append(sub_look_up_arr) 147 look_up_array = np.concatenate(look_up_arr_list)148 149 return(look_up_array)150 151def subsequence_array(look_up_array, AA_seq):152 AA_uniq = np.unique([x for x in AA_seq])153 df_list = []154 155 for AA in AA_uniq:156 sub_look_up = look_up_array[look_up_array[:,0] == AA]157 index = np.where(sub_look_up[:,3] == max(sub_look_up[:,3]))158 index = np.concatenate(index)[0]159 col_sum = sub_look_up.sum(axis = 0)[4]160 #print(sub_look_up)161 162 163 #print(len(df[:,0]))164 #print(col_sum)165 subset = []166 #print(sub_look_up)167 #print(range(len(sub_look_up))168 for i in range(len(sub_look_up)):169 subset_AA = sub_look_up[i]170 #print(subset_AA)171 if subset_AA[5] != 0:172 sequence = np.arange(start = 0, stop = subset_AA[5])173 for x in sequence:174 subset.append(subset_AA[2])175 matches = [i for i,x in enumerate(AA_seq) if x == AA]176 data = {'Index': matches,177 'Codon': subset} 178 df = pd.DataFrame(data)179 df = np.asarray(df)180 df_list.append(df)181 return(df_list)182 183def mean_dist_fitness_func_arr(array):184 distance_list = []185 for i in array:186 a = i187 array = array[1:]188 for j in array:189 if a[1] == j[1]:190 distance = j[0] - a[0]191 distance_list.append(distance)192 if len(distance_list) == 0:193 mean_distance = None194 else:195 mean_distance = statistics.mean(distance_list)196 fitness = mean_distance197 return(fitness)198 199def mutation_func(offspring, r_mut):200 if rand() < r_mut:201 index_1 = random.randint(0, offspring.shape[0] - 1)202 index_2 = random.randint(0, offspring.shape[0] - 1)203 codon_1 = offspring[index_1,1]204 codon_2 = offspring[index_2,1]205 offspring[index_1,1] = codon_2206 offspring[index_2,1] = codon_1207 208 return offspring209 210def create_initial_pop(array, n_mut, n_pop = 10):211 pop = []212 for j in range(n_pop):213 for i in range(n_mut):214 array = mutation_func(array, r_mut = 1)215 pop.append(array)216 return(pop)217 218def selection(pop, scores, k=3):219 # first random selection220 selection_ix = randint(len(pop))221 for ix in randint(0, len(pop), k-1):222 # check if better (e.g. perform a tournament)223 if scores[ix] or scores[selection_ix] is None:224 break225 if scores[ix] < scores[selection_ix]:226 selection_ix = ix227 return pop[selection_ix]228 229def genetic_algorithm(input_array, fitness_function, n_iter, n_pop, r_mut, max_s = 20):230 pop = create_initial_pop(array = input_array, n_mut = len(input_array), n_pop = n_pop)231 t_0 = time.perf_counter()232 best = pop[0]233 best_eval = fitness_function(pop[0])234 for gen in range(n_iter):235 # evaluate all candidates in the population236 scores = [fitness_function(c) for c in pop]237 # check for new best solution238 for i in range(n_pop):239 if best_eval is not None:240 if scores[i] > best_eval:241 best, best_eval = pop[i], scores[i]242 print(">New best, Generation: %s with mean: %.3f" % (gen, scores[i]))243 selected = [selection(pop, scores) for _ in range(n_pop)]244 # create the next generation245 children = list()246 for i in range(0, n_pop, 2):247 # get selected parents in pairs248 p1, p2 = selected[i], selected[i+1]249 # crossover and mutation250 for c in [p1, p2]:251 # mutation252 mutation_func(c, r_mut)253 # store for next generation254 children.append(c)255 # replace population256 pop = children257 t_n = time.perf_counter()258 #if (t_0 - t_n) > max_s:259 #break260 #print((t_n - t_0))261 return [best, best_eval]262 263def codon_opt(AA_seq, profile, n_pop = 100, n_iter = 100, fitness_function = mean_dist_fitness_func_arr):264 look_up_array = AA_look_up(AA_seq = AA_seq, profile = profile)265 sub_arrays = subsequence_array(AA_seq = AA_seq, look_up_array = look_up_array)266 complete_seq = []267 for sub_arr in sub_arrays:268 sub_seq = genetic_algorithm(input_array = sub_arr,269 fitness_function = fitness_function, 270 n_iter = n_iter, 271 n_pop = n_pop, 272 r_mut = (1.0 / float(len(sub_arr[:,0]))) * 0.5)273 complete_seq.append(sub_seq[0])274 final = np.concatenate(complete_seq, axis = 0)275 final = final[final[:, 0].argsort()]276 final = pd.DataFrame(final).rename(columns={0: "Index", 1: "Codon"})277 final = final.merge(profile.iloc[:,0:2], how = 'left', on = 'Codon')278 return(final)