Source code for ecoli.analysis.antibiotics_colony.plot
"""
Script to make majority of figures in paper.
Before running, extract colony_data.zip in the "data" folder.
Mapping between filenames in "colony_data/sim_dfs" folder and
simulation conditions in ecoli/analysis/antibiotics_colony/__init__.py.
The __init__ file also constructs a dictionary of all columns in the saved
CSV files (and gives their paths in the simulation store hierarchy).
Refer to these other scripts for the code to create remaining figures:
- ecoli/analysis/antibiotics_colony/subgen_gene_plots/
- count_subgen.py: counts for Fig. 1B
- Run with "--data data/colony_data/glc_10000_expressome.csv"
- make_fig_1b.py: create Fig. 1B
- ecoli/analysis/antibiotics_colony/snapshot_and_hist_plot.py: Fig. 2E-F
- Run with "--local data/colony_data/2022-12-08_00-35-28_562633+0000.csv"
- ecoli/analysis/antibiotics_colony/tet_dry_mass.py: Fig. 3E
- ecoli/analysis/antibiotics_colony/amp_plots.py: Fig. 4D-J, L
- Run with "--glc_data data/colony_data/2022-12-08_00-33-56_581605+0000.csv"
and "--amp_data data/colony_data/2022-12-08_17-03-56_357734+0000.csv"
- ecoli/analysis/antibiotics_colony/spatial_autocorrelation.py: Fig. S4
- Run with "data/colony_data/sim_dfs/2022-12-08_00-35-28_562633+0000.csv"
- Repeat with CSV files for other two baseline glucose simulations
to get all Moran's I and p-values in Table S1
- ecoli/analysis/proteinCountsValidation.py: Fig. S2A
- Run with "--avg_data data/colony_data/glc_10000_proteome_avgs.csv"
- ecoli/analysis/centralCarbonMetabolismScatter.py: Fig. S2B
- Run with "--numpy_data data/colony_data/glc_10000_fluxome.csv" and
"--sim_df data/colony_data/2022-12-08_00-35-28_562633+0000.csv"
"""
import ast
import argparse
import os
from itertools import combinations
import json
import matplotlib
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy.constants import N_A
from scipy.stats import spearmanr
import seaborn as sns
from typing import Any
from vivarium.library.dict_utils import deep_merge
from ecoli.analysis.antibiotics_colony import (
COUNTS_PER_FL_TO_NANOMOLAR,
DE_GENES,
EXPERIMENT_ID_MAPPING,
MAX_TIME,
SPLIT_TIME,
restrict_data,
)
from ecoli.analysis.antibiotics_colony.exploration import (
plot_exp_growth_rate,
plot_ampc_phylo,
)
from ecoli.analysis.antibiotics_colony.timeseries import (
make_tag_video,
plot_field_snapshots,
plot_tag_snapshots,
plot_timeseries,
)
from ecoli.analysis.antibiotics_colony.validation import (
plot_colony_growth,
plot_mrna_fc,
plot_protein_synth_inhib,
plot_death_timescale_analysis,
)
[docs]
def make_figure_1a(data, metadata):
# Generational (ompF) vs sub-generational (marR) expression
columns_to_plot = {
"ompF mRNA": "0.4",
"marR mRNA": (0, 0.4, 1),
"OmpF monomer": "0.4",
"MarR monomer": (0, 0.4, 1),
}
fig, axes = plt.subplots(2, 2, sharex="col", figsize=(6, 6))
axes = np.ravel(axes)
# Arbitrarily pick a surviving agent to plot trace of
highlight_agent = "011001001"
print(f"Highlighted agent: {highlight_agent}")
plot_timeseries(
data=data,
axes=axes,
columns_to_plot=columns_to_plot,
highlight_lineage=highlight_agent,
background_lineages=False,
)
axes[0].set_xlabel(None)
axes[1].set_xlabel(None)
# Put gene name on top and remove superfluous axes labels
gene_1 = axes[0].get_ylabel().split(" ")[0]
gene_2 = axes[1].get_ylabel().split(" ")[0]
axes[0].set_ylabel("mRNA\n(counts)")
axes[2].set_ylabel("Monomer\n(counts)")
axes[0].set_title(f"Exponential: {gene_1}", fontsize=10)
axes[1].set_title(f"Sub-generational: {gene_2}", fontsize=10)
axes[0].yaxis.set_label_coords(-0.3, 0.5)
axes[2].yaxis.set_label_coords(-0.3, 0.5)
axes[1].yaxis.label.set_visible(False)
axes[3].yaxis.label.set_visible(False)
axes[0].xaxis.set_visible(False)
axes[0].spines.bottom.set_visible(False)
axes[1].xaxis.set_visible(False)
axes[1].spines.bottom.set_visible(False)
for ax in axes:
[item.set_fontsize(8) for item in ax.get_xticklabels()]
[item.set_fontsize(8) for item in ax.get_yticklabels()]
ax.xaxis.label.set_fontsize(9)
ax.yaxis.label.set_fontsize(9)
ax.xaxis.set_label_coords(0.5, -0.2)
ax.tick_params(axis="both", which="major")
fig.set_size_inches(4, 3)
plt.draw()
plt.tight_layout()
plt.subplots_adjust(wspace=0.2, hspace=0.3)
os.makedirs("out/analysis/paper_figures/", exist_ok=True)
plt.savefig(
f"out/analysis/paper_figures/fig_1a_{highlight_agent}.svg", bbox_inches="tight"
)
plt.close()
print("Done with Figure 1A.")
[docs]
def make_figure_s3a(data, metadata):
# Doubling time histogram vs literature
doubling_times = []
grouped_data = data.groupby(["Seed"])
for condition, condition_data in grouped_data:
agent_ids = condition_data.loc[:, "Agent ID"].unique()
grouped_agents = condition_data.groupby(["Agent ID"])
for agent, agent_data in grouped_agents:
# Exclude cells alive at final time point because they
# have not completed their cell cycle. Also exclude cells
# that died prematurely from cell wall defects
if agent + "0" in agent_ids:
doubling_time = (
agent_data.loc[:, "Time"].max() - agent_data.loc[:, "Time"].min()
)
doubling_times.append(doubling_time)
fig, ax = plt.subplots(figsize=(2, 2))
# Literature doubling time = 44 minutes (10.1128/jb.119.1.270-281.1974)
ax.vlines(44, 0, 150, colors=["tab:orange"])
ax.text(
45,
150,
"experiment\n\u03c4 = 44 min.",
verticalalignment="top",
horizontalalignment="left",
c="tab:orange",
fontsize=8,
)
doubling_times = np.array(doubling_times) / 60
ax.hist(doubling_times)
ax.set_xlabel("Doubling time (min)")
ax.set_ylabel("# of simulated cells")
ax.set_ylim()
sim_avg = doubling_times.mean()
ax.vlines(sim_avg, ax.get_ylim()[0], 100, linestyles=["dashed"], colors=["k"])
ax.text(
sim_avg + 2,
100,
f"simulation\n\u03c4 = {np.round(sim_avg, 1)} min.",
verticalalignment="top",
horizontalalignment="left",
c="tab:blue",
fontsize=8,
)
ax.spines.right.set_visible(False)
ax.spines.top.set_visible(False)
plt.tight_layout()
os.makedirs("out/analysis/paper_figures/", exist_ok=True)
plt.savefig(
"out/analysis/paper_figures/fig_s3a_doubling_time.svg", bbox_inches="tight"
)
plt.close()
print("Done with Figure S3A.")
[docs]
def make_figure_2a(data, metadata):
# Snapshots of colony at 1.8 hour increments with environmental
# glucose concentration colored and a single lineage highlighted
# in blue
final_timestep = data.loc[data.loc[:, "Time"] == MAX_TIME, :]
agent_ids = final_timestep.loc[:, "Agent ID"]
# Arbitrarily select a lineage to highlight
highlight_agent = agent_ids[100]
print(f"Highlighted agent: {highlight_agent}")
plot_field_snapshots(
data=data,
metadata=metadata,
highlight_lineage=highlight_agent,
highlight_color=(0, 0.4, 1),
min_pct=0.8,
colorbar_decimals=2,
)
print("Done with Figure 2A.")
[docs]
def make_figure_2b_d(data, metadata):
# Timeseries of dry mass and mRNA/protein concentrations for
# OmpF, TolC, AmpC, and MarR with highlighted lineage in blue
# Use same highlighted agent as in Figure 2A
final_timestep = data.loc[data.loc[:, "Time"] == MAX_TIME, :]
agent_ids = final_timestep.loc[:, "Agent ID"]
highlight_agent = agent_ids[100]
print(f"Highlighted agent: {highlight_agent}")
# Set up subplot layout for timeseries plots
fig = plt.figure()
gs = fig.add_gridspec(3, 4)
axes = [fig.add_subplot(gs[0, :])]
for i in range(4):
axes.append(fig.add_subplot(gs[2, i]))
for i in range(4):
axes.append(fig.add_subplot(gs[1, i], sharex=axes[i + 1]))
columns_to_plot = {
"Dry mass": (0, 0.4, 1),
}
plot_timeseries(
data=data,
axes=axes,
columns_to_plot=columns_to_plot,
highlight_lineage=highlight_agent,
)
columns_to_plot = {
"OmpF monomer": (0, 0.4, 1),
"TolC monomer": (0, 0.4, 1),
"AmpC monomer": (0, 0.4, 1),
"MarR monomer": (0, 0.4, 1),
"ompF mRNA": (0, 0.4, 1),
"tolC mRNA": (0, 0.4, 1),
"ampC mRNA": (0, 0.4, 1),
"marR mRNA": (0, 0.4, 1),
}
# Use periplasmic volume for proteins that localize to periplasm
# or outer membrane
periplasmic = ["OmpF monomer", "AmpC monomer", "TolC monomer"]
for column in columns_to_plot:
if column in periplasmic:
data.loc[:, column] /= data.loc[:, "Volume"] * 0.2
else:
data.loc[:, column] /= data.loc[:, "Volume"] * 0.8
data.loc[:, column] *= COUNTS_PER_FL_TO_NANOMOLAR
monomer_cols = list(columns_to_plot)[:4]
# Convert monomer concentrations to uM
data.loc[:, monomer_cols] /= 1000
plot_timeseries(
data=data,
axes=axes[1:],
columns_to_plot=columns_to_plot,
highlight_lineage=highlight_agent,
)
# Add more regularly spaced tick marks to dry mass timeseries
time_ticks = axes[0].get_xticks()
new_ticks = np.arange(1, np.ceil(time_ticks[1]), 1).astype(int)
# No need for tick at 7 since final tick is 7.2
new_ticks = new_ticks[new_ticks != 7].tolist()
time_ticks = [0] + new_ticks + [time_ticks[1]]
axes[0].set_xticks(ticks=time_ticks, labels=time_ticks)
# Put gene name on top and remove superfluous axes labels
gene = axes[1].get_ylabel().split(" ")[0]
axes[0].set_ylabel("Dry mass (fg)")
axes[5].set_title(gene, fontsize=12, fontweight="bold")
axes[5].set_ylabel("mRNA (nM)")
axes[1].set_ylabel("Protein (\u03bcM)")
for i in range(2, 5):
gene = axes[i].get_ylabel().split(" ")[0]
axes[i].yaxis.label.set_visible(False)
axes[4 + i].set_title(gene, fontsize=12, fontweight="bold")
axes[4 + i].yaxis.label.set_visible(False)
for ax in axes[5:]:
ax.xaxis.set_visible(False)
ax.spines.bottom.set_visible(False)
for ax in axes:
[item.set_fontsize(8) for item in ax.get_xticklabels()]
[item.set_fontsize(8) for item in ax.get_yticklabels()]
ax.xaxis.label.set_fontsize(10)
ax.yaxis.label.set_fontsize(10)
ax.tick_params(axis="both", which="major")
fig.set_size_inches(7, 4)
plt.tight_layout()
plt.subplots_adjust(hspace=0.2, wspace=0.3)
for ax in axes[1:]:
ax.xaxis.set_label_coords(0.5, -0.3)
left, bottom, width, height = ax.get_position().bounds
ax.set_position((left, bottom - 0.15, width, height))
left, bottom, width, height = axes[0].get_position().bounds
axes[0].set_position((left, bottom + 0.03, width, height))
axes[0].xaxis.set_label_coords(0.5, -0.3)
axes[0].yaxis.set_label_coords(-0.09, 0.5)
axes[5].yaxis.set_label_coords(-0.5, 0.5)
axes[1].yaxis.set_label_coords(-0.5, 0.5)
# Prettify axes (moving axis titles in to save space)
for ax in axes[1:5]:
xmin, xmax = ax.get_xlim()
ax.set_xticks([(xmin + xmax) / 2], labels=[ax.get_xlabel()], minor=True)
ax.set_xlabel(None)
ax.tick_params(
which="minor",
width=0,
length=ax.xaxis.get_major_ticks()[0].get_tick_padding(),
labelsize=10,
)
# Calculate average monomer concentrations per cell then
# average across all cell
average_agent_data = data.groupby("Agent ID").mean().mean()
# Use average simulated volume to calculate literature concentrations
volume = average_agent_data["Volume"]
# 10.1016/j.cell.2014.02.033: absolute monomer counts
# for MG1655 grown in MOPS minimal media + glucose
lit_protein_concs = {
"OmpF monomer": 71798,
"TolC monomer": 3141,
"AmpC monomer": 132,
"MarR monomer": 20,
}
orange = (1, 100 / 255, 0)
for ax, monomer in zip(axes[1:5], lit_protein_concs):
monomer_vol = volume * 0.2 if monomer in periplasmic else volume * 0.8
monomer_conc = (
lit_protein_concs[monomer] / monomer_vol * COUNTS_PER_FL_TO_NANOMOLAR / 1000
)
print(monomer, monomer_conc, monomer_vol)
ax.scatter(
MAX_TIME / 3600 + 0.3, monomer_conc, c=orange, marker="_", clip_on=False
)
ax.scatter(
MAX_TIME / 3600 + 0.3,
average_agent_data[monomer],
c=(0, 0.4, 1),
marker="_",
clip_on=False,
)
ax.scatter(
MAX_TIME / 3600 + 0.5,
monomer_conc,
c=orange,
marker="o",
clip_on=False,
s=20,
)
ax.scatter(
MAX_TIME / 3600 + 0.5,
average_agent_data[monomer],
c=(0, 0.4, 1),
marker="o",
clip_on=False,
s=20,
)
# 10.1126/science.abk2066: mRNA fractions (concentration over
# total mRNA concentration) from K-12 strain NCM3722 in glucose
# minimal medium (2 biological replicates)
lit_mrna_concs = {
"ompF mRNA": (8.57e-03, 8.08e-03),
"tolC mRNA": (3.21e-04, 3.06e-04),
"ampC mRNA": (6.55e-06, 7.09e-06),
"marR mRNA": (7.88e-06, 7.13e-06),
}
mrna_counts = json.load(open("data/colony_data/glc_10000_total_mrna.json", "r"))
# Exclude cells from end of simulation and cells that died because they
# did not complete their normal cell cycle (inaccurate average mRNA count)
exclude_cells = data.loc[data["Time"] == MAX_TIME, "Agent ID"].tolist()
cracked_cells = data.loc[data["Wall cracked"], "Agent ID"].unique()
unique_agents = data["Agent ID"].unique().tolist()
for cracked_agent in cracked_cells:
if cracked_agent + "0" not in unique_agents:
exclude_cells.append(cracked_agent)
complete_mrna_counts = [
mrna_count
for agent, mrna_count in mrna_counts.items()
if agent not in exclude_cells
]
average_mrna_count = np.mean(complete_mrna_counts)
for ax, mrna in zip(axes[5:], lit_mrna_concs):
mrna_vol = volume * 0.8
mrna_conc = (
np.mean(lit_mrna_concs[mrna])
* average_mrna_count
/ mrna_vol
* COUNTS_PER_FL_TO_NANOMOLAR
)
print(mrna, mrna_conc, mrna_vol)
ax.scatter(
MAX_TIME / 3600 + 0.3, mrna_conc, c=orange, marker="_", clip_on=False
)
ax.scatter(
MAX_TIME / 3600 + 0.3,
average_agent_data[mrna],
c=(0, 0.4, 1),
marker="_",
clip_on=False,
)
ax.scatter(
MAX_TIME / 3600 + 0.5, mrna_conc, c=orange, marker="o", clip_on=False, s=20
)
ax.scatter(
MAX_TIME / 3600 + 0.5,
average_agent_data[mrna],
c=(0, 0.4, 1),
marker="o",
clip_on=False,
s=20,
)
# Bring back disappearing y ticks
for ax in axes:
yticks = ax.get_yticks()
ax.set_yticks(yticks, yticks.astype(int))
# Get fractional upper limit for MarR monomer expression
marR_max = np.round(data.loc[:, "MarR monomer"].max(), 1)
axes[4].set_ylim([0, marR_max])
axes[4].set_yticks([0, marR_max], [0, marR_max], fontsize=8)
axes[4].spines.left.set_bounds([0, marR_max])
os.makedirs("out/analysis/paper_figures/", exist_ok=True)
plt.savefig(
"out/analysis/paper_figures/fig_2b_d_timeseries.svg", bbox_inches="tight"
)
plt.close()
print("Done with Figure 2B-D.")
[docs]
def make_figure_s1(data, metadata):
# Make plots describing glucose depletion from environment
# Convert glucose uptake to units of mmol / g DCW / hr
exchange_data = data.loc[:, ["Dry mass", "Exchanges"]]
glc_flux = exchange_data.apply(
lambda x: x["Exchanges"]["GLC[p]"]
/ N_A
* 1000
/ (x["Dry mass"] * 1e-15)
/ 2
* 3600,
axis=1,
)
data["Glucose intake"] = -glc_flux
grouped_data = data.groupby("Agent ID").mean()
print(
"Mean glucose intake (mmol/g DCW/hr): "
+ str(grouped_data.loc[:, "Glucose intake"].mean())
)
print(
"Std. dev. glucose intake (mmol/g DCW/hr): "
+ str(grouped_data.loc[:, "Glucose intake"].std())
)
fig, ax = plt.subplots(figsize=(3, 3))
sns.histplot(
grouped_data.loc[:, "Glucose intake"], color=(0, 0.4, 1), ax=ax, linewidth=1
)
ax.set_xlabel("Glucose intake\n(mmol/g DCW/hr)", fontsize=10)
ax.set_ylabel("Simulated cells", fontsize=10)
ax.set_xticks(ax.get_xticks(), ax.get_xticks().astype(int), fontsize=10)
ax.set_yticks(ax.get_yticks(), ax.get_yticks().astype(int), fontsize=10)
plt.tight_layout()
os.makedirs("out/analysis/paper_figures/", exist_ok=True)
plt.savefig(
"out/analysis/paper_figures/fig_s1a_glc_intake.svg", bbox_inches="tight"
)
plt.close()
print("Done with Figure S1A.")
# Cross-section of environmental glucose concentration
field_data = metadata["Glucose"][10000]["fields"]
xticks = np.arange(0, 50, 5)
xcoords = xticks + 2.5
sample_times = [10400, 15600, 20800, 26000]
cmap = matplotlib.colormaps["Greys"]
norm = matplotlib.colors.Normalize(vmin=0, vmax=26000)
fig, ax = plt.subplots(figsize=(3, 3))
for sample_time in sample_times:
cross_sec = np.array((field_data[sample_time]["GLC[p]"])).T[4]
color = cmap(norm(sample_time))
ax.plot(xcoords, cross_sec, c=color)
ax.scatter(xcoords, cross_sec, c=color)
ax.text(
51,
cross_sec.mean(),
f"{np.round(sample_time/3600, 1)} hr",
horizontalalignment="left",
verticalalignment="center",
c=color,
)
ax.set_xlabel("Distance from left edge\nof environment (\u03bcm)")
ax.set_ylabel("Cross-sectional glucose (mM)")
xticks = np.append(xticks, 50).astype(int)
ax.set_xticks(xticks)
plt.savefig("out/analysis/paper_figures/fig_s1b_env_cross.svg", bbox_inches="tight")
plt.close()
print("Done with Figure S1B.")
[docs]
def make_figure_s5(data, metadata):
# Perform linear regression on average per-cell protein expression for
# all possible pairs of the four antibiotic resistance genes
data = restrict_data(data)
monomers = ["OmpF monomer", "AmpC monomer", "TolC monomer", "MarR monomer"]
data.loc[:, monomers] = (
data.loc[:, monomers].divide(data.loc[:, "Volume"], axis=0)
* COUNTS_PER_FL_TO_NANOMOLAR
/ 1000
)
avg_concs = data.loc[:, monomers + ["Agent ID"]].groupby("Agent ID").mean()
new_colnames = {col: col + " (\u03bcM)" for col in monomers}
avg_concs.rename(columns=new_colnames, inplace=True)
g = sns.pairplot(avg_concs, kind="reg", corner=True)
combos = list(combinations(new_colnames.values(), 2))
plot_xvars = np.array(g.x_vars)
plot_yvars = np.array(g.y_vars)
for monomer_1, monomer_2 in combos:
# Add Spearman R and p-value to corner of each plot
r, p = spearmanr(avg_concs[monomer_1], avg_concs[monomer_2])
adj_p = p * len(combos)
adj_p = min(adj_p, 1)
x_idx = np.where(plot_xvars == monomer_1)[0][0]
y_idx = np.where(plot_yvars == monomer_2)[0][0]
if x_idx > y_idx:
x_idx = np.where(plot_xvars == monomer_2)[0][0]
y_idx = np.where(plot_yvars == monomer_1)[0][0]
ax = g.axes[y_idx, x_idx]
ax.text(
s=f"r = {np.round(r, 2)}",
y=1,
x=1,
ha="right",
va="top",
transform=ax.transAxes,
)
if adj_p < 0.05:
ax.text(
s=f"p = {np.format_float_scientific(adj_p, 1)}*",
y=0.9,
x=1,
ha="right",
va="top",
transform=ax.transAxes,
weight="bold",
)
else:
ax.text(
s=f"p = {np.format_float_scientific(adj_p, 1)}",
y=0.9,
x=1,
ha="right",
va="top",
transform=ax.transAxes,
)
print(
f"{monomer_1} vs. {monomer_2}: r = {r},"
f" Bonferroni corrected p = {adj_p}"
)
os.makedirs("out/analysis/paper_figures/", exist_ok=True)
plt.savefig("out/analysis/paper_figures/fig_s5_supp_protein_pairplot.svg")
plt.close()
print("Done with Figure S5.")
[docs]
def make_figure_s3b(data, metadata):
# Plot std. dev. in average birth time of cells for each generation
agent_data = data.groupby("Agent ID")
start_times = []
for agent_id, agent in agent_data:
# Exclude final generation of cells because not all cells in the prior
# generation had finished dividing
if len(agent_id) > 9:
continue
start_times.append((len(agent_id), agent["Time"].min()))
start_times = pd.DataFrame(start_times)
start_times.rename(columns={0: "Generation", 1: "Start time (s)"}, inplace=True)
std_start_times = start_times.groupby("Generation").std()
fig, ax = plt.subplots(figsize=(2, 2))
ax.scatter(std_start_times.index, std_start_times["Start time (s)"], c="k")
ax.set_xlabel("Generation", fontsize=9)
ax.set_ylabel("Birth time std. dev. (s)", fontsize=9)
ax.set_xticks(range(2, 10), range(2, 10), fontsize=8)
ax.set_yticks(range(0, 600, 100), range(0, 600, 100), fontsize=8)
sns.despine(ax=ax, trim=True, offset=3)
ax.set_xticks(range(2, 10), range(2, 10), fontsize=8)
ax.set_yticks(range(0, 600, 100), range(0, 600, 100), fontsize=8)
os.makedirs("out/analysis/paper_figures/", exist_ok=True)
plt.savefig(
"out/analysis/paper_figures/fig_s3b_birth_time_std_dev.svg", bbox_inches="tight"
)
plt.close()
print("Done with Figure S3B.")
[docs]
def make_figure_3l(data, metadata):
# Plot total colony mass traces for all tested tetracycline concs.
fig, ax = plt.subplots(figsize=(2.1, 2.35))
plot_colony_growth(data, ax)
offset_time = SPLIT_TIME / 3600
min_time = np.round(-offset_time, 1)
max_time = np.round(MAX_TIME / 3600 - offset_time, 1)
ticks = [min_time, 0, int(max_time)]
ax.set_xticks(np.array(ticks) + offset_time, ticks, size=8)
ax.spines["bottom"].set_bounds(0, MAX_TIME / 3600)
ax.spines["left"].set_bounds(ax.get_ylim())
ax.set_xlabel("Hours after tet. addition", size=9)
ax.set_ylabel("Colony mass (fg)", size=9)
yticklabels = [f"$10^{int(exp)}$" for exp in np.log10(ax.get_yticks())]
ax.set_yticks(ax.get_yticks(), yticklabels, size=9)
plt.tight_layout()
os.makedirs("out/analysis/paper_figures/", exist_ok=True)
fig.savefig(
"out/analysis/paper_figures/fig_3l_tet_colony_mass.svg", bbox_inches="tight"
)
plt.close()
print("Done with Figure 3L.")
[docs]
def make_figure_3d_and_3m(data, metadata):
# Plot snapshots with intervals of 1.3 hours of decrease in
# instantaneous growth rate after tetracycline addition (3D)
# Also, create two scatterplots of average doubling rate
# against active ribosome concentration for cells in the first
# and fourth hours of tetracycline exposure (3M)
data = data.loc[data.loc[:, "Time"] <= MAX_TIME, :]
data = data.sort_values(["Condition", "Agent ID", "Time"])
# Draw blue border around highlighted agent lineage
highlight_agent_id = "0011111"
plot_exp_growth_rate(data, metadata, highlight_agent_id)
plt.close()
print("Done with Figure 3D.")
[docs]
def make_figure_3f_h(data, metadata):
# Timeseries traces of tetracycline and active ribosome
# concentrations in minutes surrounding tet. addition
# Filter data to only include 150 seconds before and after
glucose_mask = (
(data.loc[:, "Time"] >= 11400)
& (data.loc[:, "Time"] <= SPLIT_TIME)
& (data.loc[:, "Condition"] == "Glucose")
)
tet_mask = (
(data.loc[:, "Time"] >= SPLIT_TIME)
& (data.loc[:, "Time"] <= 11700)
& (data.loc[:, "Condition"] == "Tetracycline (1.5 mg/L)")
)
transition_data = data.loc[glucose_mask | tet_mask, :]
# Convert tetracycline concentrations to uM
transition_data.loc[:, "Periplasmic tetracycline"] *= 1000
transition_data.loc[:, "Cytoplasmic tetracycline"] *= 1000
# Convert to concentration using cytoplasmic volume
transition_data.loc[:, "Active ribosomes"] /= transition_data.loc[:, "Volume"] * 0.8
# Convert all concentrations to uM
transition_data.loc[:, "Active ribosomes"] *= COUNTS_PER_FL_TO_NANOMOLAR / 1000
transition_data.rename(
columns={
"Periplasmic tetracycline": "Tetracycline\n(periplasm)",
"Cytoplasmic tetracycline": "Tetracycline\n(cytoplasm)",
"Active ribosomes": "Active\nribosomes",
},
inplace=True,
)
fig, axes = plt.subplots(1, 3, figsize=(5, 1.5))
short_term_columns = {
"Tetracycline\n(periplasm)": 0,
"Tetracycline\n(cytoplasm)": 1,
"Active\nribosomes": 2,
}
for column, ax_idx in short_term_columns.items():
plot_timeseries(
data=transition_data,
axes=[axes.flat[ax_idx]],
columns_to_plot={column: (0, 0.4, 1)},
highlight_lineage="0011111",
filter_time=False,
background_alpha=0.5,
background_linewidth=0.3,
)
for ax in axes.flat:
ylim = ax.get_ylim()
yticks = np.round(ylim, 0).astype(int)
ax.set_yticks(yticks, yticks, size=9)
ax.set_xlabel(None)
# Mark minutes since tetracycline addition
ax.set_xticks(
ticks=[
11430 / 3600,
11490 / 3600,
11550 / 3600,
11610 / 3600,
11670 / 3600,
],
labels=[-2, -1, 0, 1, 2],
size=9,
)
ax.spines.bottom.set(
bounds=(11400 / 3600, 11700 / 3600),
linewidth=1,
visible=True,
color=(0, 0, 0),
alpha=1,
)
ylabel = ax.get_ylabel()
ax.set_ylabel(None)
ax.set_title(ylabel, size=9)
axes.flat[0].set_ylabel("\u03bcM", size=9, labelpad=0)
axes.flat[1].set_xlabel("Minutes after tetracycline addition", size=9)
# Ensure that active ribosomes plot starts at y = 0
axes.flat[-1].set_ylim(0, axes.flat[-1].get_ylim()[1])
new_yticks = [0, axes.flat[-1].get_yticks()[-1]]
axes.flat[-1].set_yticks(new_yticks, new_yticks)
axes.flat[-1].spines["left"].set_bounds(new_yticks)
plt.tight_layout()
fig.subplots_adjust(left=0.08, right=0.99, top=0.8, bottom=0.3, wspace=0.35)
os.makedirs("out/analysis/paper_figures/", exist_ok=True)
fig.savefig("out/analysis/paper_figures/fig_3f_h_tet_short_term.svg")
plt.close()
print("Done with Figure 3F-H.")
[docs]
def make_figure_3i_k(data, metadata):
# Timeseries traces of micF-ompF duplex, ompF mRNA, and OmpF monomer
# concentrations in hours surrounding tetracycline addition
# Filter data to include glucose for first 11550 seconds and
# tetracycline data for remainder of simulation
long_transition_data = restrict_data(data)
long_term_columns = {
"micF-ompF duplex": 0,
"ompF mRNA": 1,
"OmpF monomer": 2,
}
# Use periplasmic volume for OmpF monomer because it localizes
# to outer membrane
periplasmic = ["OmpF monomer"]
for column in long_term_columns:
if column in periplasmic:
long_transition_data.loc[:, column] /= (
long_transition_data.loc[:, "Volume"] * 0.2
)
else:
long_transition_data.loc[:, column] /= (
long_transition_data.loc[:, "Volume"] * 0.8
)
# Convert all concentrations to uM
long_transition_data.loc[:, column] *= COUNTS_PER_FL_TO_NANOMOLAR / 1000
fig, axes = plt.subplots(1, 3, figsize=(5, 1.5))
for column, ax_idx in long_term_columns.items():
plot_timeseries(
data=long_transition_data,
axes=[axes.flat[ax_idx]],
columns_to_plot={column: (0, 0.4, 1)},
highlight_lineage="0011111",
filter_time=False,
background_alpha=0.5,
background_linewidth=0.3,
)
split_hours = SPLIT_TIME / 3600
rounded_split_hours = np.round(split_hours, 1)
for ax in axes.flat:
ylim = ax.get_ylim()
yticks = np.round(ylim, 0).astype(int)
ax.set_yticks(yticks, yticks, size=9)
# Mark hours since tetracycline addition
xlim = np.array(ax.get_xlim())
xticks = np.append(xlim, split_hours)
xtick_labels = np.trunc(xticks - split_hours).astype(int).tolist()
xtick_labels = [
label if label != int(-split_hours) else -rounded_split_hours
for label in xtick_labels
]
ax.set_xticks(ticks=xticks, labels=xtick_labels, size=9)
ax.set_xlabel(None)
ax.spines.bottom.set(
bounds=(0, MAX_TIME / 3600),
linewidth=1,
visible=True,
color=(0, 0, 0),
alpha=1,
)
ylabel = ax.get_ylabel()
ax.set_ylabel(None)
ax.set_title(ylabel, size=9, pad=12)
axes.flat[0].set_ylabel("\u03bcM", size=9, labelpad=-6)
fig.supxlabel("Hours after tetracycline addition", size=9)
# Ensure that OmpF monomer plot starts at y = 0
for i in range(2):
y_max = np.round(axes.flat[i].get_ylim()[-1], 2)
new_yticks = [0, y_max]
axes.flat[i].set_yticks(new_yticks, new_yticks)
axes.flat[i].spines["left"].set_bounds(new_yticks)
axes.flat[-1].set_ylim(0, axes.flat[-1].get_ylim()[1])
new_yticks = [0, axes.flat[-1].get_yticks()[-1]]
axes.flat[-1].set_yticks(new_yticks, new_yticks)
axes.flat[-1].spines["left"].set_bounds(new_yticks)
plt.tight_layout()
fig.subplots_adjust(left=0.08, right=0.99, top=0.8, bottom=0.3, wspace=0.35)
os.makedirs("out/analysis/paper_figures/", exist_ok=True)
fig.savefig("out/analysis/paper_figures/fig_3i_k_tet_long_term.svg")
plt.close()
print("Done with Figure 3I-K.")
[docs]
def make_figure_s6b(data, metadata):
# Scatter plot of simulated and measured mRNA fold changes
# after exposure to 1.5 mg/L tetracycline. acrA, acrB, tolC
# genes highlighted in blue
genes_to_plot = DE_GENES.loc[:, "Gene name"]
fig, ax = plt.subplots(1, 1, figsize=(3, 3))
genes = ["acrA", "acrB", "tolC"]
highlight_genes = {gene: (0, 0.4, 1) for gene in genes}
plot_mrna_fc(data, ax, genes_to_plot, highlight_genes=highlight_genes)
ax.spines["left"].set_bounds((-1, 1.5))
ax.set_yticks([-1, -0.5, 0, 0.5, 1, 1.5])
os.makedirs("out/analysis/paper_figures/", exist_ok=True)
fig.savefig(
"out/analysis/paper_figures/fig_s6b_tet_gene_exp.svg", bbox_inches="tight"
)
plt.close()
print("Done with Figure S6B.")
[docs]
def make_figure_s6a(data, metadata):
# Dose-response curve for protein synthesis inhibition in model
# and real cells
# Jenner et al. 2013: 10.1073/pnas.1216691110
jenner = pd.read_csv("data/colony_data/jenner_2013.csv", header=None).rename(
columns={0: "Tetracycline", 1: "Percent inhibition"}
)
jenner["Source"] = ["Jenner et al. 2013"] * len(jenner)
jenner.loc[:, "Percent inhibition"] = 100 - (jenner.loc[:, "Percent inhibition"])
# Olson et al. 2006: 10.1128/AAC.01499-05
olson = pd.read_csv("data/colony_data/olson_2006.csv", header=None).rename(
columns={0: "Tetracycline", 1: "Percent inhibition"}
)
olson["Source"] = ["Olson et al. 2006"] * len(olson)
literature = pd.concat([jenner, olson])
fig, ax = plt.subplots(1, 1, figsize=(3, 3))
plot_protein_synth_inhib(data, ax, literature)
plt.tight_layout()
plt.subplots_adjust(top=1, bottom=0.25, left=0.2, right=1)
ax.set_xlabel("Tetracycline (\u03bcM)", size=10)
ax.set_ylabel("Protein synthesis inhibition (%)", size=10)
ax.tick_params(axis="both", which="major", labelsize=10)
ax.spines["bottom"].set_bounds((0.08, 300))
children = ax.get_children()
for i in [4, 5, 6]:
children[i].set_fontsize(10)
os.makedirs("out/analysis/paper_figures/", exist_ok=True)
plt.savefig(
"out/analysis/paper_figures/fig_s6a_protein_synth_inhib.svg",
bbox_inches="tight",
)
plt.close()
print("Done with Figure S6A.")
[docs]
def make_figure_s7(data, metadata):
# Decrease in outer membrane permeability but increase in
# cytoplasmic tetracycline and decreased AcrAB-TolC conc.
data = restrict_data(data)
data["Time"] -= SPLIT_TIME
column = "Outer tet. permeability (cm/s)"
# Convert permeability to nm / s for readability
data[column] *= 1e7
fig, ax = plt.subplots(figsize=(3, 3))
plot_timeseries(
data=data,
axes=[ax],
columns_to_plot={column: (0, 0.4, 1)},
highlight_lineage="0011111",
filter_time=False,
background_alpha=0.5,
background_linewidth=0.3,
)
data[column].max() * 1e7
ax.set_ylabel("OM tet. perm. (nm/s)")
ax.set_xlabel("Hours after tetracycline addition")
ax.set_xticks(np.append(ax.get_xticks(), 0), np.append(ax.get_xticks(), 0))
os.makedirs("out/analysis/paper_figures/", exist_ok=True)
plt.savefig(
"out/analysis/paper_figures/fig_s7a_tet_om_perm.svg", bbox_inches="tight"
)
plt.close()
# Get AcrAB-TolC periplasmic concentration in uM
data["AcrAB-TolC (\u03bcM)"] = (
data["AcrAB-TolC"] / (data["Volume"] * 0.2) * COUNTS_PER_FL_TO_NANOMOLAR / 1000
)
data["Time"] /= 60
# Convert tetracycline concentrations to uM
data["Initial external tet."] *= 1000
data["Cytoplasmic tetracycline"] *= 1000
data.rename(
columns={
"Cytoplasmic tetracycline": "Cytoplasmic tetracycline (\u03bcM)",
"Time": "Minutes after tetracycline addition",
"Initial external tet.": "External tet. (\u03bcM)",
},
inplace=True,
)
# Start 2 min after tet. addition to skip initial spike in tet conc.
tet_data = data.loc[
data["Minutes after tetracycline addition"] >= 2,
[
"External tet. (\u03bcM)",
"Cytoplasmic tetracycline (\u03bcM)",
"Periplasmic tetracycline (\u03bcM)",
"Minutes after tetracycline addition",
"AcrAB-TolC (\u03bcM)",
],
]
# Highlight MIC in blue and use grayscale for other concentrations
cmap = matplotlib.colormaps["Greys"]
antibiotic_min = data.loc[:, "External tet. (\u03bcM)"].min()
antibiotic_max = data.loc[:, "External tet. (\u03bcM)"].max()
norm = matplotlib.colors.Normalize(
vmin=1.5 * antibiotic_min - 0.5 * antibiotic_max, vmax=antibiotic_max
)
antibiotic_concs = data.loc[:, "External tet. (\u03bcM)"].unique()
palette = {
antibiotic_conc: cmap(norm(antibiotic_conc))
for antibiotic_conc in antibiotic_concs
}
palette[3.375] = (0, 0.4, 1)
fig, ax = plt.subplots(figsize=(3, 3))
sns.lineplot(
tet_data,
x="Minutes after tetracycline addition",
y="AcrAB-TolC (\u03bcM)",
hue="External tet. (\u03bcM)",
ax=ax,
palette=palette,
errorbar=None,
)
max_time = np.round((MAX_TIME - SPLIT_TIME) / 60, 0)
ax.set_xlim(2, int(max_time))
xticks = ax.get_xticks()[:-1]
xticks = np.insert(xticks[xticks != 0], 0, 2)
xticks = np.append(xticks, int(max_time))
ax.set_xticks(xticks, xticks.astype(int))
plt.tight_layout()
ax.legend(bbox_to_anchor=(1.1, 1), title="Ext. tet. (\u03bcM)")
sns.despine(ax=ax, trim=True, offset=3)
plt.savefig(
"out/analysis/paper_figures/fig_s7c_acrab_tet_concs.svg", bbox_inches="tight"
)
plt.close()
# Subtract avg. cytoplasmic tetracycline concentration at 2 min
# post-tetracycline addition
initial_steady_tet = (
tet_data.loc[tet_data["Minutes after tetracycline addition"] == 2, :]
.groupby("External tet. (\u03bcM)")
.mean()
)
tet_data = (
tet_data.set_index("External tet. (\u03bcM)") - initial_steady_tet
).reset_index()
fig, ax = plt.subplots(figsize=(3, 3))
sns.lineplot(
tet_data,
x="Minutes after tetracycline addition",
y="Cytoplasmic tetracycline (\u03bcM)",
ax=ax,
hue="External tet. (\u03bcM)",
palette=palette,
errorbar=None,
legend=None,
)
max_time = np.round((MAX_TIME - SPLIT_TIME) / 60, 0)
ax.set_xlim(2, int(max_time))
xticks = ax.get_xticks()[:-1]
xticks = np.insert(xticks[xticks != 0], 0, 2)
xticks = np.append(xticks, int(max_time))
ax.set_xticks(xticks, xticks.astype(int))
ax.set_ylabel("Increase in cyto. tet. (\u03bcM)")
plt.tight_layout()
sns.despine(ax=ax, trim=True, offset=3)
plt.savefig(
"out/analysis/paper_figures/fig_s7b_tet_conc_cyto.svg", bbox_inches="tight"
)
plt.close()
[docs]
def make_figure_s8(data, metadata):
# Long-term positive feedback in growth inhibition
# Filter data to include glucose for first 11550 seconds and
# tetracycline data for remainder of simulation
long_transition_data = restrict_data(data)
long_term_columns = {"Active RNAP": (0, 0.4, 1), "mRNA mass": (0, 0.4, 1)}
# Convert RNAP count to uM using cytoplasmic volume
long_transition_data.loc[:, "Active RNAP"] /= (
long_transition_data.loc[:, "Volume"] * 0.8
)
long_transition_data.loc[:, "Active RNAP"] *= COUNTS_PER_FL_TO_NANOMOLAR / 1000
fig, axes = plt.subplots(1, 2, figsize=(5, 1.5))
plot_timeseries(
data=long_transition_data,
axes=axes.flat,
columns_to_plot=long_term_columns,
highlight_lineage="0011111",
filter_time=False,
background_alpha=0.5,
background_linewidth=0.3,
)
split_hours = SPLIT_TIME / 3600
rounded_split_hours = np.round(split_hours, 1)
for ax in axes.flat:
ylim = ax.get_ylim()
yticks = np.round(ylim, 0).astype(int)
ax.set_yticks(yticks, yticks, size=9)
# Mark hours since tetracycline addition
xlim = np.array(ax.get_xlim())
xticks = np.append(xlim, split_hours)
xtick_labels = np.trunc(xticks - split_hours).astype(int).tolist()
xtick_labels = [
label if label != int(-split_hours) else -rounded_split_hours
for label in xtick_labels
]
ax.set_xticks(ticks=xticks, labels=xtick_labels, size=9)
ax.set_xlabel(None)
ax.spines.bottom.set(
bounds=(0, MAX_TIME / 3600),
linewidth=1,
visible=True,
color=(0, 0, 0),
alpha=1,
)
ylabel = ax.get_ylabel()
ax.set_ylabel(None)
ax.set_title(ylabel, size=9, pad=12)
axes.flat[0].set_ylabel("\u03bcM", size=9, labelpad=-6)
fig.supxlabel("Hours after tetracycline addition", size=9)
# Ensure that plots start at y = 0
for i in range(2):
y_max = np.round(axes.flat[i].get_ylim()[-1], 2)
new_yticks = [0, y_max]
axes.flat[i].set_yticks(new_yticks, new_yticks)
axes.flat[i].spines["left"].set_bounds(new_yticks)
axes.flat[-1].set_ylim(0, axes.flat[-1].get_ylim()[1])
new_yticks = [0, axes.flat[-1].get_yticks()[-1]]
axes.flat[-1].set_yticks(new_yticks, new_yticks)
axes.flat[-1].spines["left"].set_bounds(new_yticks)
# mRNA mass has units fg
axes[-1].set_ylabel("fg", size=9, labelpad=-6)
plt.tight_layout()
fig.subplots_adjust(left=0.08, right=0.99, top=0.8, bottom=0.3, wspace=0.35)
os.makedirs("out/analysis/paper_figures/", exist_ok=True)
fig.savefig("out/analysis/paper_figures/fig_s8_tet_supp_feedback.svg")
plt.close()
print("Done with Figure S8.")
[docs]
def make_figure_4l(data, metadata):
# Timeseries traces of total colony mass for all tested ampicillin concs.
fig, ax = plt.subplots(figsize=(2.1, 2.1))
plot_colony_growth(
data, ax, antibiotic_col="Initial external amp.", mic=5.724, antibiotic="Amp."
)
offset_time = SPLIT_TIME / 3600
min_time = np.round(-offset_time, 1)
max_time = np.round(MAX_TIME / 3600 - offset_time, 1)
ticks = [min_time, 0, int(max_time)]
ax.set_xticks(np.array(ticks) + offset_time, ticks, size=8)
ax.spines["bottom"].set_bounds(0, MAX_TIME / 3600)
ax.spines["left"].set_bounds(ax.get_ylim())
ax.set_xlabel("Hours after amp. addition", size=9)
ax.set_ylabel("Colony mass (fg)", size=9)
yticklabels = [f"$10^{int(exp)}$" for exp in np.log10(ax.get_yticks())]
ax.set_yticks(ax.get_yticks(), yticklabels, size=9)
plt.tight_layout()
os.makedirs("out/analysis/paper_figures/", exist_ok=True)
fig.savefig("out/analysis/paper_figures/fig_4l_amp_colony_mass.svg")
plt.close()
print("Done with Figure 4L.")
[docs]
def make_figure_4c(data, metadata):
# Snapshots of colony in 1.3 hour intervals after ampicillin addition
# with ratio of current gap area over average untreated gap area colored
# as log 2 fold change
# Get fold change over average glucose porosity
data["Relative porosity"] = (
data.loc[:, "Porosity"] * data.loc[:, "Extension factor"]
)
mean_glc_porosity = data.loc[
data.loc[:, "Condition"] == "Glucose", "Relative porosity"
].mean()
fc_col = "Total defect area\n($\\mathregular{log_2}$ fold change)"
data.loc[:, fc_col] = np.log2(data.loc[:, "Relative porosity"] / mean_glc_porosity)
# For 10 seconds after division, cells are assigned a porosity of 0
# until the next update from the cell wall process. Clean up these values
# so colormap works properly.
data.loc[data.loc[:, fc_col] == -np.inf, fc_col] = 0
# Set up custom divergent colormap
cmp = matplotlib.colors.LinearSegmentedColormap.from_list(
"divergent", [(0, 0, 0), (1, 1, 1), (0.678, 0, 0.125)]
)
magnitude = data.loc[:, fc_col].abs().max()
norm = matplotlib.colors.Normalize(vmin=-magnitude, vmax=magnitude)
snapshot_times = np.array([1.9, 3.2, 4.5, 5.8, 7.1]) * 3600
snapshot_times = np.array([3.2, 4.5, 5.8, 7.1]) * 3600
# Draw blue border around highlighted agent lineage
highlight_agent_id = "001111111"
highlight_agent_ids = [
highlight_agent_id[: i + 1] for i in range(len(highlight_agent_id))
]
highlight_agent = {
agent_id: {"membrane_width": 0.5, "membrane_color": (0, 0.4, 1)}
for agent_id in highlight_agent_ids
}
fig = plot_tag_snapshots(
data=data,
metadata=metadata,
tag_colors={fc_col: {"cmp": cmp, "norm": norm}},
snapshot_times=snapshot_times,
return_fig=True,
figsize=(6, 1.5),
highlight_agent=highlight_agent,
show_membrane=True,
)
fig.axes[0].set_xticklabels(
np.abs(np.round(fig.axes[0].get_xticks() / 3600 - SPLIT_TIME / 3600, 1))
)
fig.axes[0].set_xlabel("Hours after ampicillin addition")
os.makedirs("out/analysis/paper_figures/", exist_ok=True)
fig.savefig(
"out/analysis/paper_figures/fig_4c_amp_snapshots.svg", bbox_inches="tight"
)
plt.close()
print("Done with Figure 4C.")
[docs]
def make_figure_4m(data, metadata):
# Plot phylogenetic tree of cells in a representative ampicillin sim,
# where each node is a cell colored by average AmpC concentration
os.makedirs("out/analysis/paper_figures/", exist_ok=True)
plot_ampc_phylo(data)
print("Done with Figure 4M.")
[docs]
def make_figure_4n(data, metadata):
# Compare average generations to lysis in our simulations
# against literature data
fig, ax = plt.subplots(1, 1, figsize=(2.2, 2.25))
axs = [ax]
plot_death_timescale_analysis(data, axs)
os.makedirs("out/analysis/paper_figures/", exist_ok=True)
plt.savefig("out/analysis/paper_figures/fig_4n_misc.svg")
plt.close()
print("Done with Figure 4N.")
[docs]
def make_figure_s10(data, metadata):
# Pairplot of avg. AmpC, OmpF, PBP1B, and PBP1B concentration against
# avg. periplasmic ampicillin concentration (with Spearman R and p val.)
data = restrict_data(data)
data = data.loc[data.loc[:, "Time"] >= SPLIT_TIME, :]
protein_names = [
"AmpC monomer",
"OmpF monomer",
"PBP1a complex",
"PBP1b gamma complex",
"AcrAB-TolC",
]
# Convert everything into uM (use periplasmic concentration for all
# because of outer membrane/periplasmic localization)
data[protein_names] = (
data[protein_names].divide(data["Volume"] * 0.2, axis=0)
* COUNTS_PER_FL_TO_NANOMOLAR
/ 1000
)
data["Periplasmic ampicillin"] *= 1000
agent_ids = data.loc[:, "Agent ID"].unique()
final_agents = data.loc[data.loc[:, "Time"] == MAX_TIME, "Agent ID"].unique()
# Return True if cell died
dead_dict = {}
def check_death(agent_id):
status = dead_dict.get(agent_id)
if status:
return status
if agent_id + "0" not in agent_ids:
if agent_id not in final_agents:
dead_dict[agent_id] = True
return True
dead_dict[agent_id] = False
return False
data["Dead"] = data["Agent ID"].apply(check_death)
add_vars = ["Agent ID", "Dead", "Periplasmic ampicillin"]
avg_data = data[protein_names + add_vars].groupby("Agent ID").mean()
avg_data["Dead"] = avg_data["Dead"].astype(bool)
rename_cols = {curr_name: curr_name + " (\u03bcM)" for curr_name in protein_names}
rename_cols["Periplasmic ampicillin"] = "Periplasmic ampicillin (\u03bcM)"
avg_data = avg_data.rename(columns=rename_cols)
fig, axs = plt.subplots(2, 3, sharey=True)
for i, column in enumerate(protein_names):
column += " (\u03bcM)"
plot_legend = False
if i == len(protein_names) - 1:
plot_legend = True
sns.scatterplot(
avg_data,
x=column,
y="Periplasmic ampicillin (\u03bcM)",
hue="Dead",
ax=axs.flat[i],
legend=plot_legend,
)
r, p = spearmanr(avg_data[column], avg_data["Periplasmic ampicillin (\u03bcM)"])
# Bonferroni correction
p = p * len(protein_names)
if p > 1:
p = 1
print(f"{column}: r = {r}, p = {p}")
sns.despine(ax=axs.flat[i])
axs.flat[i].text(
s=f"r = {np.round(r, 2)}",
y=1,
x=1,
ha="right",
va="top",
transform=axs.flat[i].transAxes,
)
if p < 0.05:
axs.flat[i].text(
s=f"p = {np.format_float_scientific(p, 1)}*",
y=0.9,
x=1,
ha="right",
va="top",
transform=axs.flat[i].transAxes,
weight="bold",
)
else:
axs.flat[i].text(
s=f"p = {np.format_float_scientific(p, 1)}",
y=0.9,
x=1,
ha="right",
va="top",
transform=axs.flat[i].transAxes,
)
plt.tight_layout()
sns.move_legend(axs.flat[-2], "center right", bbox_to_anchor=(2.0, 0.5))
axs.flat[-1].remove()
plt.savefig("out/analysis/paper_figures/fig_s10_protein_amp_corr.svg")
plt.close()
print("Done with Figure S10.")
[docs]
def make_figure_videos(data: pd.DataFrame, metadata: dict[str, Any]):
"""
One-off video of sub-generational ampC expression and cell death at MIC.
"""
data["AmpC monomer (\u03bcM)"] = (
data["AmpC monomer"]
/ (data.loc[:, "Volume"] * 0.2)
* COUNTS_PER_FL_TO_NANOMOLAR
/ 1000
)
make_tag_video(
data=data,
metadata=metadata,
tag_colors={"AmpC monomer (\u03bcM)": (0.6, 1, 1)},
out_prefix="ampicillin_2mg_L",
)
[docs]
def load_exp_data(experiment_ids):
data = []
metadata = {}
for exp_id in experiment_ids:
exp_data = pd.read_csv(
f"data/colony_data/sim_dfs/{exp_id}.csv",
dtype={"Agent ID": str, "Seed": int},
index_col=0,
)
if exp_data.loc[:, "Dry mass"].iloc[-1] == 0:
exp_data = exp_data.iloc[:-1, :]
data.append(exp_data)
with open(f"data/colony_data/sim_dfs/{exp_id}_metadata.json", "r") as f:
metadata = deep_merge(metadata, json.load(f))
data = pd.concat(data)
# Convert strings to dictionaries
data["Boundary"] = data["Boundary"].apply(ast.literal_eval)
initial_external_tet = []
initial_external_amp = []
for condition in data["Condition"].unique():
cond_data = data.loc[data.loc[:, "Condition"] == condition, :]
if condition == "Glucose":
initial_external_tet += [0] * len(cond_data)
initial_external_amp += [0] * len(cond_data)
continue
curr_len = len(initial_external_tet)
for boundary_data in cond_data.loc[:, "Boundary"]:
# Assumes only one antibiotic is used at a time
tet_conc = boundary_data["external"]["tetracycline"]
if tet_conc != 0:
initial_external_tet += [tet_conc] * len(cond_data)
initial_external_amp += [0] * len(cond_data)
break
amp_conc = boundary_data["external"]["ampicillin[p]"]
if amp_conc != 0:
initial_external_amp += [amp_conc] * len(cond_data)
initial_external_tet += [0] * len(cond_data)
break
if len(initial_external_tet) == curr_len:
initial_external_tet += [0] * len(cond_data)
initial_external_amp += [0] * len(cond_data)
data["Initial external tet."] = initial_external_tet
data["Initial external amp."] = initial_external_amp
return data, metadata
[docs]
def main():
tet_local = [
"2022-12-08_01-13-41_036971+0000",
"2022-12-08_01-37-02_043920+0000",
"2022-12-08_01-37-17_383563+0000",
"2022-12-08_01-37-25_382616+0000",
"2022-12-08_01-37-31_999399+0000",
"2022-12-08_01-37-38_566402+0000",
"2022-12-08_01-37-44_216110+0000",
"2022-12-08_01-37-52_725211+0000",
"2022-12-08_01-37-57_809101+0000",
"2022-12-08_01-38-03_635076+0000",
"2022-12-08_01-38-09_020029+0000",
]
conditions = {
"1a": ["Glucose"],
"s3a": ["Glucose"],
"2a": ["Glucose"],
"2b_d": ["Glucose"],
"s1": ["Glucose"],
"s5": ["Glucose"],
"s3b": ["Glucose"],
"3l": [
"Glucose",
"Tetracycline (0.5 mg/L)",
"Tetracycline (1 mg/L)",
"Tetracycline (1.5 mg/L)",
"Tetracycline (2 mg/L)",
"Tetracycline (4 mg/L)",
],
"3d_and_3m": [
"Glucose",
"Tetracycline (0.5 mg/L)",
"Tetracycline (1 mg/L)",
"Tetracycline (1.5 mg/L)",
"Tetracycline (2 mg/L)",
"Tetracycline (4 mg/L)",
],
"3f_h": ["Glucose", "Tetracycline (1.5 mg/L)"],
"3i_k": ["Glucose", "Tetracycline (1.5 mg/L)"],
"s6b": ["Glucose", "Tetracycline (1.5 mg/L)"],
"s6a": [str(i) for i in range(11)],
"s7": [
"Glucose",
"Tetracycline (0.5 mg/L)",
"Tetracycline (1 mg/L)",
"Tetracycline (1.5 mg/L)",
"Tetracycline (2 mg/L)",
"Tetracycline (4 mg/L)",
],
"s8": ["Glucose", "Tetracycline (1.5 mg/L)"],
"4l": [
"Glucose",
"Ampicillin (0.5 mg/L)",
"Ampicillin (1 mg/L)",
"Ampicillin (1.5 mg/L)",
"Ampicillin (2 mg/L)",
"Ampicillin (4 mg/L)",
],
"4c": ["Glucose", "Ampicillin (2 mg/L)"],
"4m": ["Glucose", "Ampicillin (2 mg/L)"],
"4n": [
"Ampicillin (0.5 mg/L)",
"Ampicillin (1 mg/L)",
"Ampicillin (1.5 mg/L)",
"Ampicillin (2 mg/L)",
"Ampicillin (4 mg/L)",
],
"videos": ["Ampicillin (2 mg/L)"],
}
seeds = {
"1a": [10000],
"s3a": [10000],
"2a": [10000],
"2b_d": [10000],
"s1": [10000],
"s5": [10000],
"s3b": [10000],
"3l": [0],
"3d_and_3m": [0],
"3f_h": [0],
"3i_k": [0],
"s6b": [0, 100, 10000],
"s6a": [0],
"s7": [0],
"s8": [0],
"4l": [0],
"4c": [0],
"4m": [0],
"4n": [0],
"videos": [0],
}
parser = argparse.ArgumentParser()
parser.add_argument(
"--fig_ids",
"-f",
help="List of lowercase figure IDs to create (e.g. 1a). \
Default is all.",
nargs="+",
choices=list(seeds.keys()),
)
args = parser.parse_args()
if args.fig_ids is None:
args.fig_ids = conditions.keys() - {"3f"}
ids_to_load = []
for fig_id in args.fig_ids:
if fig_id == "3f":
ids_to_load.extend(tet_local)
continue
for condition in conditions[fig_id]:
for seed in seeds[fig_id]:
ids_to_load.append(EXPERIMENT_ID_MAPPING[condition][seed])
# De-duplicate IDs while preserving order
ids_to_load = list(dict.fromkeys(ids_to_load))
data, metadata = load_exp_data(ids_to_load)
for fig_id in args.fig_ids:
filter = np.isin(data.loc[:, "Condition"], conditions[fig_id]) & np.isin(
data.loc[:, "Seed"], seeds[fig_id]
)
fig_data = data.loc[filter, :].copy()
globals()[f"make_figure_{fig_id}"](fig_data, metadata)
if __name__ == "__main__":
main()