RT-PCR Pooling algorithm

RTPCR.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy\n",
    "import scipy.optimize\n",
    "np.random.seed(seed=42)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "SAMPLES=960 #Total number of samples to be checked\n",
    "RATE=0.01 # The precentage of positive samples out of the total number of samples\n",
    "WELLS=96 # Number of wells in the assay\n",
    "FRACT=0.5 # Fraction of samples (out of total) to be put in each well.\n",
    "\n",
    "S = np.zeros((SAMPLES,WELLS)) # This is the S matrix\n",
    "#Randomly choose positive samples\n",
    "positives = np.random.choice(np.arange(SAMPLES),size=int(RATE*SAMPLES),replace=False)\n",
    "\n",
    "#This is the ground truth:\n",
    "X = np.zeros(SAMPLES)\n",
    "X[positives] = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[836 477 350 893 923 261 215 334  86]\n"
     ]
    }
   ],
   "source": [
    "#The index of the positive samples:\n",
    "print(positives)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Assign samples to wells:\n",
    "arr = []\n",
    "for i in range(WELLS):\n",
    "    arr.append(np.random.choice(np.arange(SAMPLES),replace=False,size=int(FRACT*SAMPLES)))\n",
    "    S[arr[-1],i] = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Compute the signal we would have measured in the RT-PCR reaction by taking the sum of all positive samples in each well\n",
    "W = X[arr].sum(axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Solve the linear system of equations (S*x = W) using least-squares\n",
    "x = scipy.optimize.lsq_linear(S.T,W,bounds=(0,1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 86 215 261 334 350 477 836 893 923]\n"
     ]
    }
   ],
   "source": [
    "#Get the positive predictions:\n",
    "\n",
    "p_pred = np.where(np.round(x.x)>0)[0]\n",
    "print(p_pred)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 86 215 261 334 350 477 836 893 923]\n"
     ]
    }
   ],
   "source": [
    "p_true = np.sort(positives)\n",
    "print(p_true)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Success Rate: 1.000000\n"
     ]
    }
   ],
   "source": [
    "acc = p_pred == p_true\n",
    "print(\"Success Rate: %f\" % acc.mean())"
   ]
  }
 ],
 "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.8.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}