forked from igemsoftware2021/CarpinchoToeholdDesigner
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathmain.py
198 lines (144 loc) · 8.24 KB
/
main.py
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
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
########################### GETTING START########################################
######## Download the NUPACK package using the command:
######## python3 -m pip install -U nupack -f ~/Downloads/nupack-4.0.0.27/package
######## for more information go to: https://docs.nupack.org/start/
#################################################################################
########################## Import Libraries ##################################
from Bio.Seq import Seq #v1.79
from nupack import * #v4.0.027
import numpy as np #v1.21.1
from random import randint #-
import pandas as pd #v1.3.1
import matplotlib.pyplot as plt #v3.4.2
import forgi.visual.mplotlib as fvm #v2.0.2
import forgi #v2.0.2
import RNA #v0.7.5
import os #-
#Import secondary commands
from Input import *
from helper import *
#Create Results folder
path = (PATH + "/ResultsToehold")
os.makedirs(path, exist_ok=True)
#Design Conditions
#Experiment conditions from input file
#Design option specify any non-default job option from previous Wolfe17 paper
# (see more here: https://docs.nupack.org/design/#job-options)
#The Wobble option allows you to set whether the design will follow WATSON-CRICK base pair rules
my_options = DesignOptions(f_stop=0.01,
seed=0,
wobble_mutations= Wobble)
#Setting soft constrainst is an alternative to hard constraint
# which prohibits certain pattern violations
# (such as equal sequences), in which case the constraints are "softer".
# For more information: https://docs.nupack.org/design/#specify-soft-constraints
my_soft_const = [Pattern(Prevent)]
model1 = Model(material = Material,
celsius= T,
sodium= Sodium,
magnesium= Magnesium)
# Create Toehold Switch Complex
#First define the domains that are part of the toehold structure
# and then generate the re-designed toehold determined by the sequence
# and structure in TargetComplex.
Structure_def = structure_define(miRNA, Loop, Linker, Reporter, paired, unpaired)
toe_strand = TargetStrand(Structure_def[2], name='Toehold')
C1 = TargetComplex([toe_strand], Structure_def[0], name="C1")
toe_complex = complex_design(complexes = [C1],
model=model1,
soft_constraints=my_soft_const,
options=my_options)
toehold_complex = toe_complex.run(trials= Trials)
#Stop codons identification and replacing in tube design trials
results_list = []
for trials in range(Trials):
trigger = Domain(miRNA, name="trigger")
Trigger = TargetStrand([trigger], name="Trigger")
toeholdcomplex = toehold_complex[trials].to_analysis(C1)
codon = rep_stop_codons(str(toeholdcomplex), Structure_def[3],StopCodons)
print(toeholdcomplex)
if codon[1]: #replace stop codons by "VNN"
toeholdcomplex = codon[0]
#Create TargetStrand to run test tube design
# and optimizate the stop codon replaced
toehold_dom = Domain(str(toeholdcomplex), name = "toehold_dom")
toehold_strand = TargetStrand([toehold_dom], name = "toehold_strand")
# Trials for test tube design
C2 = TargetComplex([Trigger, toehold_strand], Structure_def[1], name = "C2")
t_toe = TargetTube(on_targets={C2: 1e-06},
off_targets=SetSpec(max_size=2),
name='t_toe')
my_tube_design = tube_design(tubes=[t_toe],
soft_constraints=my_soft_const,
model=model1, options=my_options)
#if have replaced stop codon, run more 10 times
tubetrials = 1
if codon[1]:
tubetrials = Trials_stop
results= my_tube_design.run(trials=tubetrials)
#Define file name and directory
path_trial = (path + "/%s") % (trials)
os.makedirs(path_trial, exist_ok=True)
for tubetri in range(tubetrials):
switchfinal = results[tubetri].to_analysis(toehold_strand)
print(".")
#Re-evaluates the presence of stop codons
if rep_stop_codons(str(switchfinal), Structure_def[3], StopCodons)[1]:
print("Tube Design Version Failed - Stop Codon detected")
continue
binding = results[tubetri].to_analysis(C2)
triggerfinal = results[tubetri].to_analysis(Trigger)
#Get MFE values for binding trigger-toehold, toehold, trigger and RBS-linker
my_mfe_trigswitch = mfe(strands=binding, model=model1)
my_mfe_switch = mfe(strands = switchfinal, model = model1)
my_mfe_trigger = mfe(strands=triggerfinal, model=model1)
my_mfe_rbslinker = mfe(strands=str(switchfinal)[paired+unpaired+8:], model=model1)
#my_complex_ensemble_defect = defect(strands=[switchfinal], structure= Structure_def[0], model = model1)
#Difference between Binding MFE and the sum of unbound trigger and toehold energies
diff_mfe = my_mfe_trigswitch[0].energy - (my_mfe_trigger[0].energy + my_mfe_switch[0].energy)
#if not "."*(unpaired+3) in str(my_mfe_switch[0].structure)[:unpaired+3]:
# print("-")
# continue
CodonOpt = "%d.%d" % (trials,tubetri)
sampled_structures = sample(strands=[str(switchfinal)], num_sample=1, model=model1)
print(sampled_structures)
ensemble_defect_binding = defect(strands=[miRNA,str(switchfinal)[:51]], structure='(21+.3)21.3(5.11)5.3', model=model1)
ensemble_defect_complex = defect(strands= [str(switchfinal)], structure=str(my_mfe_switch[0].structure), model=model1)
ensemble_defect_toeholddom = defect(strands=[str(switchfinal)[3:13]], structure = '.10', model=model1)
ensemble_defect_trigger = defect(strands= [str(miRNA)], structure=".21", model=model1)
print("complex")
print(ensemble_defect_complex)
print("binding")
print(ensemble_defect_binding)
print("single strand complex")
print(ensemble_defect_toeholddom)
print("trigger")
print(ensemble_defect_trigger)
#score = 5*ensemble_defect_trigger + 4*ensemble_defect_toeholddom + 3*ensemble_defect_complex + 1*ensemble_defect_binding
path_versions = (path_trial + "/%s") % (CodonOpt)
os.makedirs(path_versions, exist_ok=True)
probability_matrix_toehold = pairs(strands = str(switchfinal),model = model1)
probability_matrix_binding= pairs(strands=str(binding), model=model1)
save_file(CodonOpt, path_versions, str(binding), str(my_mfe_trigswitch[0].structure),
probability_matrix_binding, Structure_def[1], "Binding_Toehold_Trigger", model1)
save_file(CodonOpt, path_versions, str(switchfinal), str(my_mfe_switch[0].structure),
probability_matrix_toehold, Structure_def[0], "Toehold", model1)
#Append results
appends = [CodonOpt, miRNA,str(Seq(miRNA).reverse_complement()), str(switchfinal),
str(my_mfe_switch[0].structure), str(my_mfe_trigger[0].structure),
str(my_mfe_trigswitch[0].structure),str(my_mfe_trigger[0].energy),
str(my_mfe_switch[0].energy), str(my_mfe_trigswitch[0].energy),
str(my_mfe_rbslinker[0].energy), str(diff_mfe), ensemble_defect_complex,
ensemble_defect_binding #, CodonOpt
]
results_list.append(appends)
results_tab = pd.DataFrame(results_list, columns=["Trial", "miRNA.sequence",
"miRNA.Toehold.Complement", "Full.Toehold.seq",
"Toehold.structure", "Trigger.structure",
"Full.TriggerToehold.structure",
"Trigger.energy", "Toehold.energy",
"BindingToeholdTrigger.energy", "DeltaG_RBSLinker",
"Diff_energy", "ensemble_defect_complex", "ensemble_defect_binding"]).set_index("Trial")
print(results_tab)
results_tab.to_csv(r'Toehold_output_RBS_Bsubtilis.csv',index = True, header = True)
print("Saved and Finished")