fig2abcd_mapping_production_fitness_landscape.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# Mapping the 7-mer production fitness landscape"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# Setup"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "ExecuteTime": {
          "end_time": "2021-06-15T17:13:58.134959Z",
          "start_time": "2021-06-15T17:13:57.550003Z"
        }
      },
      "outputs": [],
      "source": [
        "import pandas as pd\n",
        "import numpy as np\n",
        "import matplotlib as mpl\n",
        "import matplotlib.pyplot as plt\n",
        "import seaborn as sns\n",
        "\n",
        "from IPython.display import Image\n",
        "\n",
        "from scipy.stats import gaussian_kde, norm\n",
        "\n",
        "from utils_f4f import heatmap, si_format"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "ExecuteTime": {
          "end_time": "2021-06-15T17:13:58.138092Z",
          "start_time": "2021-06-15T17:13:58.136004Z"
        }
      },
      "outputs": [],
      "source": [
        "mpl.rcParams['svg.fonttype'] = 'none'\n",
        "mpl.rcParams['font.family'] = ['sans-serif']\n",
        "mpl.rcParams['font.sans-serif'] = ['Arial']\n",
        "mpl.rcParams['text.usetex'] = False\n",
        "mpl.rcParams['pdf.fonttype'] = 42\n",
        "mpl.rcParams['ps.fonttype'] = 42"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "-------\n",
        "\n",
        "# Codon usage does not affect production fitness \n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "ExecuteTime": {
          "end_time": "2021-06-15T17:59:36.070118Z",
          "start_time": "2021-06-15T17:59:35.876289Z"
        }
      },
      "outputs": [],
      "source": [
        "# Data\n",
        "\n",
        "# Load production fitness scores for the two codon replicates\n",
        "df2a = pd.read_csv('data/modeling_library_production_fitness.csv')\n",
        "df2a.set_index('AA', inplace=True)\n",
        "#df2a = df2a.loc[df2a['Label'] == 'Designed']"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "ExecuteTime": {
          "end_time": "2021-06-15T17:59:36.846461Z",
          "start_time": "2021-06-15T17:59:36.830554Z"
        }
      },
      "outputs": [],
      "source": [
        "# Preprocessing \n",
        "\n",
        "# Remove non-detected in DNA (NA)\n",
        "df2a = df2a.loc[(~df2a['Production1'].isna()) & (~df2a['Production2'].isna()), :]\n",
        "\n",
        "# log2 of the production (virus to DNA enrichment)\n",
        "x = np.log2(df2a['Production1'])\n",
        "y = np.log2(df2a['Production2'])\n",
        "\n",
        "# Variants not detected in both codon replicates \n",
        "remove = np.isinf(x) & np.isinf(y)\n",
        "x = x[~remove]\n",
        "y = y[~remove]\n",
        "\n",
        "# Variants detected in both codon replicates\n",
        "keep = (~np.isinf(x)) & (~np.isinf(y))\n",
        "x_both = x[keep]\n",
        "y_both = y[keep]\n",
        "\n",
        "# Orphan codon replicates (detected in one codon replicate but not the other)\n",
        "y_missing = y[np.isinf(x)]\n",
        "x_missing = x[np.isinf(y)]\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "# Figure \n",
        "\n",
        "# Prepare plotting kernal \n",
        "kernel = gaussian_kde(np.vstack([\n",
        "    x_both.sample(n=10000, random_state=1), \n",
        "    y_both.sample(n=10000, random_state=1)\n",
        "]))\n",
        "c_both = kernel(np.vstack([x_both.values, y_both.values]))\n",
        "\n",
        "# Figure specifications \n",
        "fig = plt.figure(figsize=(1.4,1.3), dpi=150)\n",
        "\n",
        "gs = fig.add_gridspec(2, 2, left=0.275, right=0.95, bottom=0.22, top=0.95, \n",
        "                      width_ratios=[1, 6], height_ratios=[6, 1], hspace=0., wspace=0)\n",
        "ax = fig.add_subplot(gs[0, 1])\n",
        "ax.scatter(x_both, y_both, c=c_both, cmap=mpl.cm.inferno, s=0.2, edgecolor='none', rasterized=True)\n",
        "ax.set_xticks([]); ax.set_yticks([])\n",
        "xlim = [-13, 9]\n",
        "bins = np.linspace(*xlim, 25)\n",
        "ax.set_xlim(xlim); ax.set_ylim(xlim)\n",
        "ax.text(0.03, 0.97, r'$r$ = {:.2f}'.format(np.corrcoef(x_both, y_both)[0, 1]),\n",
        "       transform=ax.transAxes, ha='left', va='top', fontsize=7)\n",
        "\n",
        "ax.text(0.97, 0.01, 'n={}'.format(si_format(len(x_both)), precision=2, format_str='{value}{prefix}',), \n",
        "        transform=ax.transAxes, ha='right', va='bottom', fontsize=7)\n",
        "\n",
        "# Left Marginal\n",
        "ax = fig.add_subplot(gs[0, 0])\n",
        "ax.hist(y_missing, bins=bins, edgecolor='none', orientation='horizontal', density=True, color='r')\n",
        "ax.set_ylim(xlim)\n",
        "ax.set_xticks([]); ax.set_yticks([-10, -5, 0, 5, 10, 15])\n",
        "ax.text(0.97, 0.97, 'n={}'.format(si_format(len(y_missing)), precision=2, format_str='{value}{prefix}',), \n",
        "        transform=ax.transAxes, ha='right', va='top', fontsize=7, rotation=90, color='r')\n",
        "ax.set_ylabel('Codon Rep 2', labelpad=1)\n",
        "ax.tick_params(axis='both', labelsize=7, length=2, pad=1)\n",
        "\n",
        "# Bottom Marginal\n",
        "ax = fig.add_subplot(gs[1, 1])\n",
        "ax.hist(x_missing, bins=bins, edgecolor='none', density=True, color='r')\n",
        "ax.set_xlim(xlim)\n",
        "ax.set_xticks([-10, -5, 0, 5, 10, 15]); ax.set_yticks([])\n",
        "ax.text(0.97, 0.8, 'n={}'.format(si_format(len(x_missing)), precision=2, format_str='{value}{prefix}',), \n",
        "        transform=ax.transAxes, ha='right', va='top', fontsize=7, color='r')\n",
        "ax.set_xlabel('Codon Rep 1', labelpad=2)\n",
        "ax.tick_params(axis='both', labelsize=7, length=2, pad=1)\n",
        "\n",
        "# Missing label (not detected in either)\n",
        "ax = fig.add_subplot(gs[1, 0])\n",
        "for spine in ['bottom', 'left']:\n",
        "    ax.spines[spine].set_visible(False)\n",
        "ax.set_xticks([]); ax.set_yticks([])\n",
        "\n",
        "ax.text(0.8, 0.8, 'Missing', transform=ax.transAxes, color='r', ha='right', va='top', fontsize=7, clip_on=False)\n",
        "\n",
        "# Save Figure and display \n",
        "filename = 'figures/fig2a_codon_usage'\n",
        "fig.savefig('{}_600dpi.svg'.format(filename), dpi=600)\n",
        "fig.savefig('{}_1200dpi.svg'.format(filename), dpi=1200)\n",
        "fig.savefig('{}.png'.format(filename), dpi=300)\n",
        "\n",
        "plt.close()\n",
        "\n",
        "Image(filename+'.png') "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "------\n",
        "\n",
        "\n",
        "\n",
        "# The 7mer Production fitness landscape is bimodal\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "ExecuteTime": {
          "end_time": "2021-07-01T19:14:03.998827Z",
          "start_time": "2021-07-01T19:14:01.900445Z"
        }
      },
      "outputs": [],
      "source": [
        "# Data\n",
        "\n",
        "# Load production fitness scores for the modeling library variants \n",
        "df2b = pd.read_csv('data/modeling_library_production_fitness.csv')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "ExecuteTime": {
          "end_time": "2021-04-28T02:19:58.402492Z",
          "start_time": "2021-04-28T02:19:58.346414Z"
        }
      },
      "source": [
        "\n",
        "- Low Gaussian: mean = -7.3, variance = 4.25, N = 35,776\n",
        "- High Gaussian: mean = 0.22, variance = 3.65, N = 28,400\n",
        "- Non-detected (at virus level) = 5,281 (~7.6%)\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "# Non-detected varaints (Inf)\n",
        "df2b['log2Enr'] = np.log2(df2b['Production'])\n",
        "valid_x = df2b['log2Enr'][~np.isinf(df2b['log2Enr'])]\n",
        "n_invalid_x = np.isinf(df2b['log2Enr']).sum()\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "ExecuteTime": {
          "end_time": "2021-07-01T19:21:57.832913Z",
          "start_time": "2021-07-01T19:21:42.244295Z"
        }
      },
      "outputs": [],
      "source": [
        "# Figure \n",
        "\n",
        "# Figure specifications \n",
        "bin_width = 0.15\n",
        "bins = np.arange(-200, 200, bin_width)\n",
        "\n",
        "sns.set_theme(style='ticks', font_scale=0.75, rc={\n",
        "    'svg.fonttype': 'none',\n",
        "    'font.sans-serif': ['Arial'],\n",
        "    'font.family': 'sans-serif',\n",
        "    'text.usetex': False,\n",
        "    'pdf.fonttype': 42,\n",
        "    'ps.fonttype': 42,\n",
        "    'font.size': 9,\n",
        "    'axes.labelsize': 9,\n",
        "    'axes.titlesize': 9,\n",
        "    'lines.linewidth': 0.5,\n",
        "    'axes.linewidth': 0.5,\n",
        "    'legend.fontsize': 9,\n",
        "    'legend.title_fontsize': 9,\n",
        "    'xtick.major.size': 3,\n",
        "    'xtick.major.pad': 3,\n",
        "    'xtick.major.width': 0.5,\n",
        "    'ytick.major.size': 3,\n",
        "    'ytick.major.pad': 3,\n",
        "    'ytick.major.width': 0.5,\n",
        "})\n",
        "fig = plt.figure(figsize=(1.3, 1.3), dpi=300)\n",
        "gs = fig.add_gridspec(\n",
        "    1, 1,\n",
        "    wspace=0, hspace=0,\n",
        "    left=0.12, right=0.95, bottom=0.25, top=0.95\n",
        ")\n",
        "\n",
        "\n",
        "ax = fig.add_subplot(gs[0, 0])\n",
        "ax.hist(valid_x, bins=bins, linewidth=0, density=True, color='#AAA')\n",
        "\n",
        "# Distributions\n",
        "xs = np.linspace(-14, 7, 1000)\n",
        "lo_n = 35_776\n",
        "hi_n = 28_400\n",
        "sum_n = lo_n + hi_n\n",
        "lo_y = norm.pdf(xs, loc=-7.3, scale=np.sqrt(4.25)) * (lo_n / sum_n)\n",
        "hi_y = norm.pdf(xs, loc=0.22, scale=np.sqrt(3.65)) * (hi_n / sum_n)\n",
        "\n",
        "\n",
        "ax.plot(xs, lo_y, '-r', linewidth=0.5, label='Low')\n",
        "ax.fill_between(\n",
        "    xs, \n",
        "    np.zeros(len(xs)), \n",
        "    y2=lo_y, \n",
        "    zorder=len(xs)+1, \n",
        "    color='r',\n",
        "    alpha=0.3\n",
        ")\n",
        "\n",
        "ax.plot(xs, hi_y, '-b', linewidth=0.5, label='High')\n",
        "ax.fill_between(\n",
        "    xs, \n",
        "    np.zeros(len(xs)), \n",
        "    y2=hi_y, \n",
        "    zorder=len(xs)+1, \n",
        "    color='b',\n",
        "    alpha=0.3\n",
        ")\n",
        "\n",
        "# Labeling \n",
        "ax.set_xlim([-18, 7])\n",
        "ax.set_xticks([-10, -5, 0, 5])\n",
        "ax.set_ylim([0, 0.13])\n",
        "ax.set_yticks([])\n",
        "\n",
        "ax.set_xlabel('Production Fitness', labelpad=0)\n",
        "ax.set_ylabel('Density', labelpad=3)\n",
        "\n",
        "\n",
        "ax.text(-7.3, 0.97, 'Low', transform=mpl.transforms.blended_transform_factory(ax.transData, ax.transAxes), \n",
        "        ha='center', va='top', color='r', fontsize=9)\n",
        "ax.text(0.22, 0.97, 'High', transform=mpl.transforms.blended_transform_factory(ax.transData, ax.transAxes), \n",
        "        ha='center', va='top', color='b', fontsize=9)\n",
        "\n",
        "ax.axvline(-14, linewidth=0.5, linestyle='-', color='#AAA')\n",
        "ax.text(\n",
        "    0.05, 0.5, '{:.1%} undetected'.format(n_invalid_x / len(df2b)),\n",
        "    ha='left', va='center', color='#AAA', transform=ax.transAxes,\n",
        "    rotation=90, fontsize=8\n",
        ")\n",
        "\n",
        "\n",
        "# Save Figure and display \n",
        "filename = 'figures/fig2b_production_fitness_space'\n",
        "fig.savefig('{}.png'.format(filename), dpi=300)\n",
        "fig.savefig('{}_600dpi.svg'.format(filename), dpi=600)\n",
        "fig.savefig('{}_1200dpi.svg'.format(filename), dpi=1200)\n",
        "\n",
        "plt.close()\n",
        "\n",
        "Image(filename + '.png')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "-----\n",
        "\n",
        "# NNK fit sequence space is biased compared to production-fit space \n",
        "\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "# Data\n",
        "\n",
        "# NNK \n",
        "nnk = pd.read_csv('data/nnk_library_top_production_fitness_75k.csv')\n",
        "\n",
        "# Designed \n",
        "modeling = pd.read_csv('data/modeling_library_production_fitness.csv')\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "# Preprocessing \n",
        "\n",
        "# Slicing out high production fit variants \n",
        "high_fit = modeling[(np.log2(modeling.Production) > -2)] \n",
        "high_fit = high_fit.reset_index(drop = True)\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "# Figure \n",
        "\n",
        "# Configurations\n",
        "cmap = mpl.cm.coolwarm\n",
        "vmin = -0.1\n",
        "vmax = 0.1\n",
        "\n",
        "sns.set_theme(style='white', font_scale=0.75, rc={\n",
        "    'svg.fonttype': 'none',\n",
        "    'font.sans-serif': ['Arial'],\n",
        "    'font.family': 'sans-serif',\n",
        "    'text.usetex': False,\n",
        "    'pdf.fonttype': 42,\n",
        "    'ps.fonttype': 42,\n",
        "    'font.size': 7,\n",
        "    'axes.labelsize': 7,\n",
        "    'axes.titlesize': 7,\n",
        "    'lines.linewidth': 0.5,\n",
        "    'axes.linewidth': 0.5,\n",
        "    'legend.fontsize': 9,\n",
        "    'legend.title_fontsize': 9,\n",
        "    'xtick.major.size': 3,\n",
        "    'xtick.major.pad': 3,\n",
        "    'xtick.major.width': 0.5,\n",
        "    'ytick.major.size': 3,\n",
        "    'ytick.major.pad': 3,\n",
        "    'ytick.major.width': 0.5,\n",
        "})\n",
        "\n",
        "fig = plt.figure(figsize=(1.7, 1.6), dpi=300)\n",
        "gs = fig.add_gridspec(\n",
        "    1, 2,\n",
        "    hspace=0., wspace=0.15,\n",
        "    left=0.02, right=0.71, bottom=0.05, top=0.85,\n",
        "    width_ratios=[1, 1]\n",
        ")\n",
        "gs2 = fig.add_gridspec(\n",
        "    1, 1,\n",
        "    left=0.71, right=0.745, bottom=0.05, top=0.85,\n",
        ")\n",
        "\n",
        "\n",
        "# Plot heatmaps \n",
        "ax = fig.add_subplot(gs[:, 0])\n",
        "hm1 = heatmap(nnk['AA'], ax, 'NNK High\\nN=75K','results/nnk_high_75k_heatmap.csv', vmin, cmap, vmax)\n",
        "\n",
        "ax = fig.add_subplot(gs[:, 1])\n",
        "hm2 = heatmap(high_fit['AA'], ax, 'Modeling High\\nN=26K','results/modeling_high_26k_heatmap.csv',vmin, cmap, vmax)\n",
        "\n",
        "\n",
        "# Colorbar\n",
        "ax = fig.add_subplot(gs2[0, 0])\n",
        "colornorm = mpl.colors.Normalize(vmin=vmin, vmax=vmax)\n",
        "cb1 = mpl.colorbar.ColorbarBase(\n",
        "    ax, \n",
        "    cmap=cmap, norm=colornorm,\n",
        "    values=np.linspace(-0.05, 0.1, 100),\n",
        "    orientation='vertical'\n",
        ")\n",
        "cb1.set_label(\n",
        "    'AA frequency relative to expected', \n",
        "    rotation=-90, labelpad=9, fontsize=7\n",
        ")\n",
        "ax.yaxis.set_ticks_position('right')\n",
        "ax.yaxis.set_label_position('right')\n",
        "    \n",
        "ax.tick_params(axis='y', length=2, pad=1, labelsize=7)\n",
        "\n",
        "# Save Figure and display\n",
        "filename = 'figures/fig2c_production_fitness_space_aa_distribution'\n",
        "fig.savefig('{}.png'.format(filename))\n",
        "fig.savefig('{}_600dpi.svg'.format(filename),dpi=600)\n",
        "fig.savefig('{}_1200dpi.svg'.format(filename),dpi=1200)\n",
        "\n",
        "plt.close()\n",
        "\n",
        "Image(filename + '.png')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "-----\n",
        "\n",
        "# Production fitness is replicable across libraries \n",
        "\n",
        "\n",
        "Production fitness replication quality of the control set (10K) that is shared between the modeling and assessment libraries."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "# Data \n",
        "\n",
        "# Load variants of the Modeling and Assess libraries\n",
        "modeling = pd.read_csv('data/modeling_library_production_fitness.csv')\n",
        "assess = pd.read_csv('data/assessment_library_production_fitness.csv')\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "# Preprocessing\n",
        "\n",
        "# Log2 Enrichment\n",
        "modeling['log2Enr'] = np.log2(modeling['Production'])\n",
        "assess['log2Enr'] = np.log2(assess['Production'])\n",
        "\n",
        "# Shared variants between the two libraries \n",
        "df2d = (\n",
        "    modeling\n",
        "    .set_index('AA')\n",
        "    [['log2Enr']]\n",
        "    .rename(columns={'log2Enr': 'modeling'})\n",
        "    .join(\n",
        "        assess\n",
        "        .set_index('AA')\n",
        "        [['log2Enr']]\n",
        "        .rename(columns={'log2Enr': 'assessment'}),\n",
        "        how='inner'\n",
        "    )\n",
        ")\n",
        "\n",
        "# Non-detected in DNA (NAs)\n",
        "df2d = df2d.loc[(~df2d['modeling'].isna()) & (~df2d['assessment'].isna()), :]\n",
        "\n",
        "df2d.head()\n",
        "\n",
        "\n",
        "# Non-detected in virus (Infs) \n",
        "x = df2d['modeling']\n",
        "y = df2d['assessment']\n",
        "\n",
        "remove = np.isinf(x) & np.isinf(y)\n",
        "x = x[~remove]\n",
        "y = y[~remove]\n",
        "\n",
        "keep = (~np.isinf(x)) & (~np.isinf(y))\n",
        "x_both = x[keep]\n",
        "y_both = y[keep]\n",
        "\n",
        "y_missing = y[np.isinf(x)]\n",
        "x_missing = x[np.isinf(y)]"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "ExecuteTime": {
          "end_time": "2021-06-15T18:02:16.945563Z",
          "start_time": "2021-06-15T18:02:16.941532Z"
        }
      },
      "outputs": [],
      "source": [
        "# Correlation between production fitness measurements in the two libraries \n",
        "np.corrcoef(x_both, y_both)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "# Figure\n",
        "\n",
        "# Plotting kernel\n",
        "kernel = gaussian_kde(np.vstack([x_both.values, y_both.values]))\n",
        "c_both = kernel(np.vstack([x_both.values, y_both.values]))\n",
        "\n",
        "# Figure specifications \n",
        "fig = plt.figure(figsize=(1.4,1.3), dpi=150)\n",
        "\n",
        "gs = fig.add_gridspec(2, 2, left=0.275, right=0.95, bottom=0.22, top=0.95, \n",
        "                      width_ratios=[1, 6], height_ratios=[6, 1], hspace=0., wspace=0)\n",
        "ax = fig.add_subplot(gs[0, 1])\n",
        "# ax.set_aspect('equal', 'box')\n",
        "ax.scatter(x_both, y_both, c=c_both, cmap=mpl.cm.inferno, s=0.2, edgecolor='none', rasterized=True)\n",
        "ax.set_xticks([]); ax.set_yticks([])\n",
        "xlim = [-13, 9]\n",
        "bins = np.linspace(*xlim, 25)\n",
        "ax.set_xlim(xlim); ax.set_ylim(xlim)\n",
        "# ax.set_title('Fc-Ly6A')\n",
        "ax.text(0.03, 0.97, r'$r$ = {:.2f}'.format(np.corrcoef(x_both, y_both)[0, 1]),\n",
        "       transform=ax.transAxes, ha='left', va='top', fontsize=7)\n",
        "\n",
        "ax.text(0.97, 0.01, 'n={}'.format(si_format(len(x_both)), precision=2, format_str='{value}{prefix}',), \n",
        "        transform=ax.transAxes, ha='right', va='bottom', fontsize=7)\n",
        "\n",
        "# Left Marginal\n",
        "ax = fig.add_subplot(gs[0, 0])\n",
        "ax.hist(y_missing, bins=bins, edgecolor='none', orientation='horizontal', density=True, color='r')\n",
        "ax.set_ylim(xlim)\n",
        "ax.set_xticks([]); ax.set_yticks([-10, -5, 0, 5, 10, 15])\n",
        "ax.text(0.5, 0.97, 'n={}'.format(len(y_missing)), \n",
        "        transform=ax.transAxes, ha='center', va='top', fontsize=7, rotation=90, color='r')\n",
        "ax.set_ylabel('Assessment Library', labelpad=1)\n",
        "ax.tick_params(axis='both', labelsize=7, length=2, pad=1)\n",
        "\n",
        "# Bottom Marginal\n",
        "ax = fig.add_subplot(gs[1, 1])\n",
        "ax.hist(x_missing, bins=bins, edgecolor='none', density=True, color='r')\n",
        "ax.set_xlim(xlim)\n",
        "ax.set_xticks([-10, -5, 0, 5, 10, 15]); ax.set_yticks([])\n",
        "ax.text(0.97, 0.8, 'n={}'.format(len(x_missing)), \n",
        "        transform=ax.transAxes, ha='right', va='top', fontsize=7, color='r')\n",
        "ax.set_xlabel('Modeling Library', labelpad=2)\n",
        "ax.tick_params(axis='both', labelsize=7, length=2, pad=1)\n",
        "\n",
        "# Missing label\n",
        "ax = fig.add_subplot(gs[1, 0])\n",
        "for spine in ['bottom', 'left']:\n",
        "    ax.spines[spine].set_visible(False)\n",
        "ax.set_xticks([]); ax.set_yticks([])\n",
        "\n",
        "ax.text(0.8, 0.8, 'Missing', transform=ax.transAxes, color='r', ha='right', va='top', fontsize=7, clip_on=False)\n",
        "\n",
        "# Save and display figure \n",
        "filename = 'figures/fig2d_production_fitness_replicability'\n",
        "fig.savefig('{}.png'.format(filename), dpi=300)\n",
        "fig.savefig('{}_600dpi.svg'.format(filename), dpi=600)\n",
        "fig.savefig('{}_1200dpi.svg'.format(filename), dpi=1200)\n",
        "\n",
        "plt.close()\n",
        "\n",
        "Image(filename+'.png') "
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3 (ipykernel)",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.9.17"
    },
    "toc": {
      "base_numbering": 1,
      "nav_menu": {},
      "number_sections": true,
      "sideBar": true,
      "skip_h1_title": false,
      "title_cell": "Table of Contents",
      "title_sidebar": "Contents",
      "toc_cell": true,
      "toc_position": {
        "height": "calc(100% - 180px)",
        "left": "10px",
        "top": "150px",
        "width": "320px"
      },
      "toc_section_display": true,
      "toc_window_display": true
    }
  },
  "nbformat": 4,
  "nbformat_minor": 4
}