1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142 | import pickle
import numpy as np
from keras.models import model_from_json
import csv
def loop_one_hot_encode(loop_seq, loop_struct):
'''
Encodes a ribozyme loop into a 3-dimensional array meant to represent its 3-domensional physical form. The
dimensions are as follows:
width - 2, to represent the nucleotides as the directly emerge from a paired stem
height - 15, to encapsulate long loops. Most positions will be blank for shorter loops
channels - 8, to represent the 4 nucleotides, each in a bonded or unbonded state
:param loop_seq: String of nucleotides representing the sequence of the loop
:param loop_struct: String of dot-bracket notation representing the structure of the loop
:return: numpy array of size (2, 1, 15, 8)
'''
encoded_loop = np.zeros((2, 1, 15, 8))
loop_codes = ['A', 'A', 'U', 'U', 'C', 'C', 'G', 'G']
if len(loop_seq) <= 30:
# Iterates forward from the 5' end of the loop and backwards from the 3' end simultaneously, marking at each
# position the correct channel.
for i in range(len(loop_seq) // 2):
encoded_loop[0, 0, i, loop_codes.index(loop_seq[i]) + (loop_struct[i] == '(' or loop_struct[i] == ')')] = 1
encoded_loop[1, 0, i, loop_codes.index(loop_seq[-(i + 1)]) + (loop_struct[-(i + 1)] == '(' or loop_struct[-(i + 1)] == ')')] = 1
# If an odd number of nucleotides, put the on at the apex of the loop into the 1st column.
if len(loop_seq) % 2 == 1:
i = len(loop_seq) // 2
encoded_loop[0, 0, i, loop_codes.index(loop_seq[i]) + (loop_struct[i] == '(' or loop_struct[i] == ')')] = 1
else:
# If the loop is too large to fit, just does the first 15 on each side.
for i in range(15):
encoded_loop[0, 0, i, loop_codes.index(loop_seq[i]) + (loop_struct[i] == '(' or loop_struct[i] == ')')] = 1
encoded_loop[1, 0, i, loop_codes.index(loop_seq[-(i + 1)]) + (loop_struct[-(i + 1)] == '(' or loop_struct[-(i + 1)] == ')')] = 1
return encoded_loop
def struct_dict_to_array(in_dict):
'''
Turns a dictionary of ribozyme sequences and structures into a 5-dimensional numpy array, containing 4-dimensional
arrays representing paired ribozyme loops for each sequence, along with an array containing the basal
gene-regulatory activity for each sequence.
:param in_dict: Dictionary of ribozyme sequences where the keys are tuples like so: (sequence string, dot-bracket
structure string) and the value is the basal gene-regulatory activity for that sequence.
:return: Tuple containing 3 arrays:
5-dimensional numpy array representing paired ribozyme loops for each sequence
1-dimensional numpy array containing the basal gene-regulatory activity for each sequence
1-dimensional array of tuples, containing the sequences of the 2 loops for each sequence
'''
out_X = []
out_y = []
out_loops = []
first, mid, last = 'GCUGUC', 'CUGAUGA', 'GAAACAGC'
# Iterates through each sequence in the dictionary
for seq in in_dict.keys():
test_seq = seq[0]
test_struct = seq[1]
# Gets the sequence and array of the first loop
loop_seq = test_seq[test_seq.find(first) + len(first): test_seq.find(mid)]
loop_struct = test_struct[test_seq.find(first) + len(first): test_seq.find(mid)]
l1 = loop_one_hot_encode(loop_seq, loop_struct)
l1_seq = loop_seq
# Gets the sequence and array of the second loop
loop_seq = test_seq[test_seq.find(mid) + len(mid): test_seq.find(last)]
loop_struct = test_struct[test_seq.find(mid) + len(mid): test_seq.find(last)]
l2 = loop_one_hot_encode(loop_seq, loop_struct)
# Pairs the 2 loops up and stores paired array and value
out_X.append(np.concatenate((l1, l2), axis=1))
out_y.append(in_dict[seq][0])
out_loops.append((l1_seq, loop_seq))
out_X = np.array(out_X).astype('float32')
out_y = np.expand_dims(out_y, axis=1)
return out_X, out_y, out_loops
# Reads in test data
test_list = pickle.load(open('Candidate_list_RNAs_min_structures.pkl', 'rb'))
rest_dict = {}
test_dict = {}
for seq in test_list:
target_dict = test_dict
target_dict[seq[0]] = seq
test_seq_dict = {}
# Pulls test data into dictionary for later conversion to array
for i in test_dict.keys():
test_seq_dict[(tuple(test_dict[i][1][0]), tuple(test_dict[i][1][1]), test_dict[i][0], test_dict[i][2])] = \
[1]
# Segments training data by structure, creating a dictionaries where each sequence in that dictionary has the same
# structure
test_segmented_dict = {}
for loop in test_seq_dict:
if (len(loop[0][0]) / 2, len(loop[0][1]) / 2, loop[1][0], loop[1][1]) not in test_segmented_dict:
test_segmented_dict[(int(len(loop[0][0]) / 2), int(len(loop[0][1]) / 2), loop[1][0], loop[1][1])] = {}
test_segmented_dict[(len(loop[0][0]) / 2, len(loop[0][1]) / 2, loop[1][0], loop[1][1])][(loop[2], loop[3])] = \
test_seq_dict[loop]
# List of structures in the segmented dictionary that have over 50 entries
test_populous_loops = [i for i in test_segmented_dict.keys() if len(test_segmented_dict[i]) > 1]
all_pr = []
all_loops = []
# For each structure segment, finds the appropriate model, pulls it, and gets predictions for sequences in that segment
for te_seg in test_segmented_dict:
try:
teX, teY, teloops = struct_dict_to_array(test_segmented_dict[te_seg])
json_file = open('Models/' + str(list(te_seg)) + ".json", 'r')
loaded_model_json = json_file.read()
json_file.close()
loaded_model = model_from_json(loaded_model_json)
loaded_model.load_weights('Models/' + str(list(te_seg)) + "model.h5")
pr = loaded_model.predict(teX, batch_size=32)
all_loops.extend(teloops)
all_pr.extend(pr)
print("Model for " + str(te_seg) + " found and used.")
except:
print("Model for " + str(te_seg) + " not found.")
# Writes the predicted values out to a csv in order of lowest predicted basal gene-regulatory activity to highest.
with open('predictions.csv', 'w', newline='') as csvfile:
writer = csv.writer(csvfile, delimiter=',',
quotechar='|', quoting=csv.QUOTE_MINIMAL)
writer.writerow(['Loop I seq', 'Loop II seq', 'Predicted basal log10(GFP/mCh)'])
best_pr = [i[0] for i in sorted(enumerate(all_pr), key=lambda x:x[1])]
for i in best_pr:
writer.writerow([all_loops[i][0], all_loops[i][1], all_pr[i][0]])
|