Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Makes a GIF showing DNA Chisel's algorithm in action
# Warning: quick and dirty code.
import dnachisel as dc
import proglog
proglog.notebook()
import matplotlib.pyplot as plt
from dna_features_viewer import BiopythonTranslator, GraphicFeature
counter = [0]
def feature_key(f):
return str(f.label).split(",")[0]
def plot_problem(problem):
fig, ax = plt.subplots(1, figsize=(12, 4))
default_props = dict(thickness=20, box_color=None,
fontdict=dict(size=14, color="white", weight="bold"))
rec = problem.to_record()
translator = BiopythonTranslator(features_properties=lambda f: default_props)
grec = translator.translate_record(rec)
_, (levels, _) = grec.plot(annotate_inline=True, ax=ax, with_ruler=False)
levels = {f.label: level for f, level in levels.items()}
grec.plot_sequence(ax=ax, background=None)
rec = problem.to_record()
evals = problem.constraints_evaluations()
rec.features = evals.locations_as_features(with_specifications=False)
default_props = dict(label=None, color='#ff9999', linewidth=0)
translator = BiopythonTranslator(features_properties=lambda f: default_props)
grec = translator.translate_record(rec)
grec.plot(ax=ax, with_ruler=False, level_offset=-1.5)
ax.set_ylim((-2, 3))
for f in problem.sequence_edits_as_features():
ax.fill_between([f.location.start - 0.5, f.location.end-0.5], -0.5, -.9,
facecolor="orange", alpha=0.15)
fontdict = dict(weight="bold", size=12)
ax.text(-1, 1, "Constraints", ha="right", va="center", fontdict=fontdict)
ax.text(-1, -0.7, "Sequence", ha="right", va="center", fontdict=fontdict)
ax.text(-1, -1.5, "Breaches", ha="right", va="center", fontdict=fontdict)
return ax, levels
def plot_stored(stored):
location, problem, local_pb = (stored[f] for f in ("location", "problem", "local_problem"))
local_rec = local_pb.to_record()
evals = local_pb.constraints_evaluations()
local_constraints_as_features = evals.success_and_failures_as_features()
L = len(problem.sequence)
ax, levels = plot_problem(problem)
start, end = location.start - 0.5, location.end - 0.5
fog_start = min(f.location.start for f in local_rec.features) - 0.5
fog_end = max(f.location.end for f in local_rec.features) - 0.5
ax.fill_between([-1, fog_start], -.3, 3, facecolor="white", alpha=0.8)
ax.fill_between([fog_end, L+1], -.3, 3, facecolor="white", alpha=0.8)
for x in fog_start, fog_end:
ax.axvline(x, ls=":", c="grey")
if (counter[0] + 1) % 2:
for x in start, end:
ax.axvline(x, ls="--", ymin=0.05, ymax=0.31, c="grey")
ax.fill_between([start, end], -0.5, -.9, facecolor="orange", alpha=0.5)
default_props = dict(label=None, color='#99ff99', linewidth=1)
translator = BiopythonTranslator(features_properties=lambda f: default_props)
grec = translator.translate_record(local_rec)
for f in local_constraints_as_features:
color = "#77dd77" if (f.qualifiers["passes"] == "true") else "#dd7777"
label = f.qualifiers["label"]
gf = GraphicFeature(f.location.start, f.location.end, color=color, label=None,
thickness=16)
grec.plot_feature(ax, gf, level=levels[label])
return ax
def plot_png(ax):
ax.figure.savefig("./pics/captures/cap_%04d.png" % counter[0], dpi=300)
ax.figure.savefig("./pics/captures_svg/cap_%04d.svg" % counter[0])
counter[0] += 1
plt.close(ax.figure)
class CustomBarLogger(proglog.TqdmProgressBarLogger):
def store_callback(self, **kw):
if "local_problem" in kw:
ax = plot_stored(self.stored)
plot_png(ax)
def bars_callback(self, bar, attr, value, old_value=None):
if all([bar == "location", attr=="index",
value==1 or value != self.bars[bar]["total"], self.stored]):
ax = plot_stored(self.stored)
plot_png(ax)
seq = "AATGCTATGCGTCGTTTATGCAATGACATGCATGCGTCGTTTATCGAATGCA"
problem = dc.DnaOptimizationProblem(
seq,
constraints = [dc.AvoidPattern("ATGC", location=(0, 31)),
dc.AvoidPattern("CG", location=(10, 45)),
dc.EnforceGCContent(mini=0.3, maxi=0.75, window=10, location=(15, 50)),
],
logger = CustomBarLogger()
)
ax, _ = plot_problem(problem)
plot_png(ax)
plot_png(ax)
problem.resolve_constraints()
import moviepy.editor as mpe
def fl(gf, t):
if int(t) % 2:
return gf(t)
p = (t % 1)**3
return p * gf(int(t+1)) + (1 - p) * gf(t)
clip = (mpe.ImageSequenceClip("./pics/captures/", fps=1)
.fl(fl)
.fx(mpe.vfx.freeze, t="end", freeze_duration=3)
.subclip(0 , -1))
clip.write_videofile("animation.mp4", fps=24)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment