fig4c_multifunction_hit_rates.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fig 4c - MultiFunction Hit rate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Imports \n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "from IPython.display import Image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Data\n",
    "Uniform = pd.read_csv('data/modeling_library_production_fitness.csv')\n",
    "F4F = pd.read_csv('data/fit4function_library_screens.csv')\n",
    "MuliFxn = pd.read_csv('data/multifunction_library.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Meta data \n",
    "Assays = ['Production','Liver','HepG2_bind','HepG2_tr','THLE_bind','THLE_tr']\n",
    "Threshold_percentage  = 50; \n",
    "Margin_std = 2;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#%% [1] Computing the hit rate for MultiFunction MultiFxnDesigneded library from validation assays  \n",
    "# Subsets\n",
    "PositiveControl = MuliFxn[MuliFxn['Label'] == 'Positive Control']\n",
    "MultiFxnDesigned = MuliFxn[MuliFxn['Label'] == 'MultiFunction']\n",
    "\n",
    "# Iterate through assays to exclude variants not matching the positive control distributions\n",
    "InxExclude = set()\n",
    "\n",
    "for assay in Assays:\n",
    "    crnt_PositiveControl = np.log2(PositiveControl[assay].values)\n",
    "    crnt_MultiFxnDesigned = np.log2(MultiFxnDesigned[assay].values)\n",
    "    \n",
    "    # Filtering out -inf values from crnt_PositiveControl\n",
    "    valid_indices = ~np.isinf(crnt_PositiveControl)\n",
    "    valid_crnt_PositiveControl = crnt_PositiveControl[valid_indices]\n",
    "    \n",
    "    mean_PositiveControl = np.nanmean(valid_crnt_PositiveControl)\n",
    "    std_PositiveControl = np.nanstd(valid_crnt_PositiveControl)\n",
    "    th_PositiveControl = max(mean_PositiveControl - Margin_std * std_PositiveControl, valid_crnt_PositiveControl.min())\n",
    "    \n",
    "    # Exclude those lower than the Threshold\n",
    "    exclude_indices = np.where(crnt_MultiFxnDesigned < th_PositiveControl)\n",
    "    InxExclude.update(exclude_indices[0])\n",
    "    \n",
    "    # Exclude those that are -inf\n",
    "    inf_indices = np.where(np.isinf(crnt_MultiFxnDesigned))\n",
    "    InxExclude.update(inf_indices[0])\n",
    "\n",
    "# Calculate Hit_rate_MultiFxn\n",
    "Hit_rate_MultiFxn = (len(MultiFxnDesigned) - len(InxExclude)) / len(MultiFxnDesigned)\n",
    "print(Hit_rate_MultiFxn)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#%% [2] Estimating the hit rate for Production-fit space from Fit4Function assays \n",
    "\n",
    "# Assuming F4F is a pandas DataFrame\n",
    "F4F['Production'] = F4F['Production2']\n",
    "\n",
    "Thresholds = []\n",
    "for assay in Assays:\n",
    "    crnt = F4F[assay].values\n",
    "    \n",
    "    # Filter out infinite and NaN values\n",
    "    valid_crnt = crnt[~np.isinf(crnt) & ~np.isnan(crnt)]\n",
    "    \n",
    "    Thresholds.append(np.percentile(valid_crnt, Threshold_percentage))\n",
    "\n",
    "# Computing for Fit4Function \n",
    "INX = set(range(len(F4F)))\n",
    "\n",
    "for assay, threshold in zip(Assays,Thresholds):\n",
    "    crnt = F4F[assay].values\n",
    "    valid_indices = np.where((crnt > threshold) & ~np.isinf(crnt) & ~np.isnan(crnt))\n",
    "    \n",
    "    # Intersect INX with valid_indices\n",
    "    INX = INX.intersection(valid_indices[0])\n",
    "\n",
    "Hit_rate_F4F = len(INX) / len(F4F)\n",
    "print(Hit_rate_F4F)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#%% [3] Estimating the hit rate for the uniform space from Fit4Function assays \n",
    "# Production fitness distribution \n",
    "# Uncomment if you want to see the bimodal distributions\n",
    "# X = np.log2(Uniform['Production'])\n",
    "# X = X[~np.isinf(X)]\n",
    "# plt.figure()\n",
    "# plt.hist(X, bins=100) \n",
    "# plt.xlabel('Production Fitness (log2)')\n",
    "# plt.ylabel('Frequency')\n",
    "# plt.show()\n",
    "\n",
    "# Processing and Calculations\n",
    "GMM_midpoint = -2.5\n",
    "Uniform['Production_log2'] = np.log2(Uniform['Production'])\n",
    "Uniform = Uniform[Uniform['Label'] == 'Designed']\n",
    "\n",
    "valid_values = Uniform['Production_log2'][(Uniform['Production_log2'] > GMM_midpoint) & (~np.isinf(Uniform['Production_log2']))]\n",
    "Production_fit_space = len(valid_values) / len(Uniform)\n",
    "#print(Production_fit_space)\n",
    "\n",
    "Hit_rate_Uniform = Hit_rate_F4F * Production_fit_space \n",
    "print(Hit_rate_Uniform)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Meta data\n",
    "\n",
    "df = pd.DataFrame()\n",
    "df['Uniform'] = np.array([Hit_rate_Uniform])\n",
    "df['HighProdFit'] = Hit_rate_F4F\n",
    "df['Designed'] = Hit_rate_MultiFxn\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Figure \n",
    "\n",
    "# Figure Configurations\n",
    "_new_black = '#000'\n",
    "sns.set_theme(style='ticks', font_scale=0.75, rc={\n",
    "    'font.family': 'sans-serif',\n",
    "    'font.sans-serif': ['Arial', 'DejaVu Sans'],\n",
    "    'svg.fonttype': 'none',\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",
    "    'axes.labelpad': 2,\n",
    "    'axes.linewidth': 0.5,\n",
    "    'axes.titlepad': 4,\n",
    "    'lines.linewidth': 0.5,\n",
    "    'legend.fontsize': 9,\n",
    "    'legend.title_fontsize': 9,\n",
    "    'xtick.labelsize': 9,\n",
    "    'ytick.labelsize': 9,\n",
    "    'xtick.major.size': 2,\n",
    "    'xtick.major.pad': 1,\n",
    "    'xtick.major.width': 0.5,\n",
    "    'ytick.major.size': 2,\n",
    "    'ytick.major.pad': 1,\n",
    "    'ytick.major.width': 0.5,\n",
    "    'xtick.minor.size': 2,\n",
    "    'xtick.minor.pad': 1,\n",
    "    'xtick.minor.width': 0.5,\n",
    "    'ytick.minor.size': 2,\n",
    "    'ytick.minor.pad': 1,\n",
    "    'ytick.minor.width': 0.5,\n",
    "\n",
    "    # Avoid black unless necessary\n",
    "    'text.color': _new_black,\n",
    "    'patch.edgecolor': _new_black,\n",
    "    'patch.force_edgecolor': False, # Seaborn turns on edgecolors for histograms by default and I don't like it\n",
    "    'hatch.color': _new_black,\n",
    "    'axes.edgecolor': _new_black,\n",
    "    # 'axes.titlecolor': _new_black # should fallback to text.color\n",
    "    'axes.labelcolor': _new_black,\n",
    "    'xtick.color': _new_black,\n",
    "    'ytick.color': _new_black\n",
    "\n",
    "})\n",
    "\n",
    "\n",
    "# Plot\n",
    "fig = plt.figure(figsize=(2.3, 1.5), dpi=150)\n",
    "gs = fig.add_gridspec(1, 1, left=0.2, right=0.95, top=0.9, bottom=0.25)\n",
    "ax = fig.add_subplot(gs[0, 0])\n",
    "\n",
    "y = np.arange(0, 1)\n",
    "\n",
    "bar_height = 0.3\n",
    "ax.barh(y - bar_height*0, df['Uniform'], height=bar_height, label='Control', color='#B3B3B3')\n",
    "ax.barh(y - bar_height*1, df['HighProdFit'], height=bar_height, label='Fit4Fxn', color='#FF5E66')\n",
    "ax.barh(y - bar_height*2, df['Designed'], height=bar_height, label='MultiFxn', color='#488ABA')\n",
    "\n",
    "ax.set_xticks(np.arange(0, 1, 0.2))\n",
    "ax.set_xticklabels([str(l) + '%' for l in np.arange(0, 100, 20)])\n",
    "ax.tick_params(axis='x', length=2, pad=2.5)\n",
    "\n",
    "ax.set_yticks(y-bar_height*1)\n",
    "ax.set_yticklabels(['Liver'])\n",
    "# ax.tick_params(axis='y', length=0, pad=2, rotation=90)\n",
    "ax.tick_params(axis='y', length=0, pad=4)\n",
    "for ticklabel in ax.get_yticklabels():\n",
    "    ticklabel.set_va('center')\n",
    "\n",
    "ax.grid(axis='x', color='#DDD', linewidth=0.5)\n",
    "\n",
    "# Labeling \n",
    "ax.set_xlabel('Hit Rate — All Functions', labelpad=4)\n",
    "\n",
    "label_x_offset = 0.02\n",
    "ax.text(df.at[0, 'Uniform'] + label_x_offset, 0, 'Uniform ({:.1%})'.format(df.at[0, 'Uniform']), \n",
    "        ha='left', va='center', color='#888', fontsize=9)\n",
    "ax.text(df.at[0, 'HighProdFit'] + label_x_offset, 0-bar_height, 'Fit4Function ({:.1%})'.format(df.at[0, 'HighProdFit']), \n",
    "        ha='left', va='center', color='#FF5E66', fontsize=9)\n",
    "ax.text(df.at[0, 'Designed'] - label_x_offset, 0-(bar_height*2)-0.01, 'MultiFunction ({:.1%})'.format(df.at[0, 'Designed']), \n",
    "        ha='right', va='center', color='#FFF', fontsize=9)\n",
    "\n",
    "# Save and show \n",
    "filename = 'figures/fig4c_multifunction_hit_rates'\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')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "f4f",
   "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"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}