ILP_formulation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Single oligo problem formulation\n",
    "\n",
    "We first present the simplified case, in which we can encode the entire library using only a single length-unconstrained oligo. \n",
    "To address this problem, we introduce the objective function $$\\underset{O}{\\text{maximize}}\\  \\sum_i t_i$$ where $1 \\leq i \\leq |T|$ is an indicator variable denoting whether target sequence $t_i \\in T$ is present in the set of genes $G$ created by the sequence library.\n",
    "\n",
    "We also introduce the variable $x_{jkl}$ where $1 \\leq j \\leq |O|$ denotes the index of oligo $o_j$ in the oligo pool $O$, $1 \\leq k \\leq L$ denotes the position of the redundant codon in the length $L$ of the gene sequence, and $1 \\leq l \\leq |C|$ denotes the use of redundant codon $c_l \\in C$ at the specified position.\n",
    "Note that for this example case $j = 1$ as we are only considering the special case of having a single oligo encode the entire library.\n",
    "Upon the variable $x_{jkl}$, we introduce the following constraint:\n",
    "\\begin{equation*}\n",
    "\\begin{aligned}\n",
    "    \\sum_l x_{jkl} & = 1 \\qquad & 1 \\leq j \\leq |O|,\\ 1 \\leq k \\leq L\\\\\n",
    "\\end{aligned}\n",
    "\\end{equation*}\n",
    "so that only a single redundant codon can be employed at each position on each oligo.\n",
    "\n",
    "Furthermore, we introduce the variable $z_{jkm}$ and the matrix $a_{lm}$, where  $1 \\leq m \\leq 22$ denotes whether the amino acid, deletion, or stop is encoded by the redundant codon indicated by index $l$.\n",
    "Therefore the following relationship exists between $x$ and $z$: $$z_{jkm} = \\sum_l x_{jkl}a_{lm}$$\n",
    "\n",
    "To test whether a specific target sequence $t_i$ can be encoded by the oligo pool, we introduce the variable $b_{ij}$.\n",
    "$b_{ij}$ indicates whether oligo $j$ can encode the target sequence $i$.\n",
    "Upon this variable we introduce the following constraints, where $s$ is defined such that $s_{ikm}$ denotes amino acid/deletion/stop $m$ at position $k$ of target sequence $i$:\n",
    "\n",
    "\\begin{equation*}\n",
    "\\begin{aligned}\n",
    "    \\sum_k\\sum_ms_{ikm}z_{jkm} - L + (L + 1)(1 - b_{ij}) & \\leq L & 1 \\leq i \\leq |T|,\\ 1 \\leq j \\leq |O| \\\\\n",
    "    \\sum_k\\sum_ms_{ikm}z_{jkm} - L + (L + 1)(1 - b_{ij}) & \\geq 0 & 1 \\leq i \\leq |T|,\\ 1 \\leq j \\leq |O| \\\\\n",
    "\\end{aligned}\n",
    "\\end{equation*}\n",
    "\n",
    "Finally, we can impose the following constraints to solve $t_i$ for all values of $i$:\n",
    "\n",
    "\\begin{equation*}\n",
    "\\begin{aligned}\n",
    "    - \\sum_jb_{ij} + (|O| + 1)t_i & \\leq |O| & 1 \\leq i \\leq |T| \\\\\n",
    "    - \\sum_jb_{ij} + (|O| + 1)t_i & \\geq 0 & 1 \\leq i \\leq |T| \\\\\n",
    "\\end{aligned}\n",
    "\\end{equation*}\n",
    "\n",
    "We note that the calculation of total produced library size $G$ requires multiplication of the total possible number of incorporated residues at each position. Because multiplcation of variables is a non-linear operation, we introduce the following constraint to account for the technology-imposed diversity limit (in the case of a single oligo, $j = 1$ only):\n",
    "$$$$\n",
    "\\begin{equation*}\n",
    "\\begin{aligned}\n",
    "    \\sum_k \\sum_l x_{jkl} \\log \\left(\\sum_m a_{lm}\\right) & \\leq \\log\\left(M\\right) & j = 1,\\ 1 \\leq k \\leq L, 1 \\leq l \\leq |C| \\\\\n",
    "\\end{aligned}\n",
    "\\end{equation*}\n",
    "where $M$ denotes the technology-imposed diversity limit. The operation $\\sum_k \\sum_l x_{jkl} x_{jkl} \\sum_m a_{lm}$ is ordinarily necessary to obtain the count of  unique residues at each position. However, the log operation cannot be applied to this value as it would result in a non-linear operation over the variable $x$. Instead, we use Lemma 1 to reformat this operation to $\\sum_k \\sum_l x_{jkl} \\log \\left(\\sum_m a_{lm}\\right)$.\n",
    "\n",
    "**Lemma 1:** For a variable $x$ and a constant matrix $a$, $\\log\\left(\\sum_i x_i \\sum_j a_{ij}\\right) = \\sum_i x_i \\log\\left(\\sum_j a_{ij}\\right)$ is true if $x_i \\in \\{0, 1\\}\\ \\forall i$, $\\sum_i x_i = 1$, and $\\sum_j a_{ij} \\geq 0\\ \\forall i$.\n",
    "\n",
    "**Proof:** Let $\\hat i$ indicate the index for which $x_{\\hat i} = 1$:\n",
    "\n",
    "\\begin{equation*}\n",
    "\\begin{aligned}\n",
    "    \\sum_i x_i \\log\\left(\\sum_j a_{ij}\\right) & = \\log\\left(\\sum_i x_i \\sum_j a_{ij}\\right) \\\\\n",
    "    \\sum_{i \\not = \\hat i} x_i \\log\\left(\\sum_j a_{ij}\\right) + x_{\\hat i} \\log\\left(\\sum_j a_{\\hat ij}\\right) & =\n",
    "    \\log\\left(\\sum_{i \\not = \\hat i} x_i \\sum_j a_{ij} + x_{\\hat i} \\sum_j a_{\\hat ij}\\right) \\\\\n",
    "    x_{\\hat i} \\log\\left(\\sum_j a_{\\hat ij}\\right) & =\n",
    "    \\log\\left(x_{\\hat i} \\sum_j a_{\\hat ij}\\right) \\\\\n",
    "    1\\log\\left(\\sum_j a_{\\hat ij}\\right) & =\n",
    "    \\log\\left(1\\right) + \\log\\left(\\sum_j a_{\\hat ij}\\right) \\\\\n",
    "    \\log\\left(\\sum_j a_{\\hat ij}\\right) & =\n",
    "    \\log\\left(\\sum_j a_{\\hat ij}\\right) \\\\\n",
    "\\end{aligned}\n",
    "\\end{equation*}\n",
    "\n",
    "\n",
    "Together, these variables, constants, and constraints give us the formulation for the ILP for a single oligo encoding the entire library:\n",
    "\n",
    "\\begin{equation*}\n",
    "\\begin{aligned}\n",
    "& \\underset{O}{\\text{maximize}}\n",
    "& & \\sum{t_i} \\\\\n",
    "& \\text{subject to}\n",
    "& & \\sum_l x_{jkl} = 1 & \\forall j,k \\\\\n",
    "& & & \\sum_l x_{jkl}a_{lm} = z_{jkm} & \\forall j,k,m \\\\\n",
    "& & & \\sum_k\\sum_ms_{ikm}z_{jkm} - L + (L + 1)(1 - b_{ij}) \\leq L & \\forall i,j,m \\\\\n",
    "& & &  \\sum_k\\sum_ms_{ikm}z_{jkm} - L + (L + 1)(1 - b_{ij}) \\geq 0 & \\forall i,j,m \\\\\n",
    "& & & - \\sum_j b_{ij} + (|O| + 1)t_i \\leq |O| & \\forall i \\\\\n",
    "& & & - \\sum_j b_{ij} + (|O| + 1)t_i \\geq 0 & \\forall i \\\\\n",
    "& & & \\sum_k \\sum_l x_{jkl} \\log \\left(\\sum_m a_{lm}\\right) \\leq \\log\\left(M\\right) & \\forall k,l,m\n",
    "\\end{aligned}\n",
    "\\end{equation*}\n",
    "\n",
    "# Multiple oligo extension"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Implementation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cvxpy\n",
    "import numpy as np\n",
    "import seq_utils"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Smallest test set"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sequences = ['ACGY',\n",
    "             'AAGY',\n",
    "             'AGAY',\n",
    "             'GGAT']\n",
    "\n",
    "n_targets, n_var_pos, S = seq_utils.create_S(sequences)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Random sequences with correlation structure"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sequences = seq_utils.generate_seqs(300, 10)\n",
    "n_targets, n_var_pos, S = seq_utils.create_S(sequences)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Valentina dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import csv\n",
    "\n",
    "sequences = []\n",
    "\n",
    "with open('del_olmo_toro_DBD_alignment.csv') as csvfile:\n",
    "    readCSV = csv.reader(csvfile, delimiter=',')\n",
    "    for row_num, row in enumerate(readCSV):\n",
    "        if row_num == 0:\n",
    "            wt = row[1:]\n",
    "        \n",
    "        seq = [i if i != '.' else wt[idx] for idx, i in enumerate(row[1:])]\n",
    "        sequences.append(''.join(seq))\n",
    "        \n",
    "n_targets, n_var_pos, S = seq_utils.create_S(sequences)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Run ILP"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "solution = seq_utils.solve_library(S, lib_lim=100, n_oligos=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t, Z, X, B = seq_utils.extract_solution(solution)\n",
    "codons = seq_utils.retrieve_codons(X)\n",
    "num_g = seq_utils.calc_num_seqs(codons)\n",
    "seqs_in = seq_utils.check_if_seqs_in_soln(sequences, codons)\n",
    "p_on_target = seq_utils.calc_prob_on_target(sequences, codons)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "codons"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_g"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sum(seqs_in)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p_on_target"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## SwiftLib comparison"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "i=0\n",
    "for seq in sequences:\n",
    "    print('>Seq{}'.format(i))\n",
    "    print(seq)\n",
    "    i+=1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sl_codons = ['ASA', 'ATA', 'NNS', 'ATA', 'WAC', 'NNS', 'GCA', 'AGM']\n",
    "sl_num_g = seq_utils.calc_num_seqs(codons)\n",
    "sl_seqs_in = seq_utils.check_if_seqs_in_soln(sequences, codons)\n",
    "sl_p_on_target = seq_utils.calc_prob_on_target(sequences, codons)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sl_num_g"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sum(sl_seqs_in)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sl_p_on_target"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}