Monte-Carlo_Simulation_for_BFG-PCA.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ideal-ancient",
   "metadata": {},
   "outputs": [],
   "source": [
    "# coding: utf-8\n",
    "import ast\n",
    "import os\n",
    "import sys\n",
    "import re\n",
    "from pprint import pprint\n",
    "import numpy as np\n",
    "from tqdm import tqdm\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "\n",
    "from func import simulation as sim\n",
    "\n",
    "##############################################\n",
    "# Monte-Carlo simulation of BFG-Y2H screen   #\n",
    "# Initially developped by Nozomu Yachie      #\n",
    "# Last modified by Daniel Evans-Yamamoto     #\n",
    "##############################################\n",
    "\n",
    "run_name = 'test'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "blessed-letter",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define teh number of strains here;\n",
    "x_strain_num = 500\n",
    "y_strain_num = 500"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "premier-discount",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Setting haploid population\n",
    "x_hap = [[\"x_%4d\"%(i+1),1] for i in range(x_strain_num)]\n",
    "y_hap = [[\"y_%4d\"%(i+1),1] for i in range(y_strain_num)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "super-tribe",
   "metadata": {},
   "outputs": [],
   "source": [
    "hap_size = 10 * 3*10**7; # 1 OD unit of haploid cells = 3 x 10^7 \n",
    "\n",
    "\n",
    "# Adjusting the cell number of each strain in pool.\n",
    "x_hap = sim.adjust_num(x_hap,0.3,hap_size)\n",
    "y_hap = sim.adjust_num(y_hap,0.3,hap_size)\n",
    "    \n",
    "cov = sim.hap_coverage(x_hap,y_hap,run_name)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "confirmed-barrier",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"The ratio of strains having at lesat the number of cells above threshold.\")\n",
    "print(\"Cells per haploid pool: %d OD600nm units\"%(hap_size/(3*10**7)))\n",
    "sim.print_LL(cov)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "scheduled-trace",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Simulating diploid cells after en masse yeast mating subjected for diploid selection\n",
    "\n",
    "mating_efficacy   = 0.01 # This is underestimation to be safe.\n",
    "diploid_selection = 500 * 3 * 10**7 # Amount of cells to be selected for diploids.\n",
    "est_diploids      = diploid_selection*mating_efficacy\n",
    "\n",
    "dip = sim.mating(x_hap,y_hap)\n",
    "dip = sim.adjust_num(dip,0.5,est_diploids)#CV of X-Y pair-dependent mating efficiencies (log-normal distribution) is 50%\n",
    "cov = sim.dip_coverage(dip,run_name)\n",
    "\n",
    "\n",
    "print(\"The ratio of diploid strains having at least the number of cells above threshold.\")\n",
    "sim.print_LL(cov)\n",
    "\n",
    "#If you want to estimate the optimal paramater for coverage, create a for loop with the variable, and plot a line plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "amateur-cabinet",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "competent-stress",
   "metadata": {},
   "outputs": [],
   "source": [
    "###############################################################\n",
    "# Diploid selection and amount of cells spread per plate\n",
    "\n",
    "\n",
    "# Amount of cells to be spread to each selection plate\n",
    "dip_plate = 50 * 1.0 * 10**7 #1 OD(600 nm) unit for diploid yeast is 1 × 10^7 cells ml^-1\n",
    "###############################################################################################################\n",
    "### !!!! After observing CFUs on the plates, update dip_plate to check if the CFU is enough down stream !!  ###\n",
    "###############################################################################################################\n",
    "\n",
    "dip = sim.adjust_num(dip,0.5,dip_plate) #CV of X-Y pair-dependent growth amplitudes (log-normal distribution) in liquid media is 50%\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eight-federation",
   "metadata": {},
   "outputs": [],
   "source": [
    "#pprint(dip)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "weird-brooks",
   "metadata": {},
   "outputs": [],
   "source": [
    "cells_for_dox_induction = 25 * 1.0 * 10**7\n",
    "dip_neg = sim.adjust_num(dip,0.1,cells_for_dox_induction) #CV of X-Y pair-dependent growth amplitudes (log-normal distribution) in non-selectable condition is 10%\n",
    "\n",
    "pos_rate  = 0.001 #PPI ositive rate is 0.1% (overestimation)\n",
    "CV_autoactivity = 0.1 # Log CV for autoactivity in BFG-PCA is 10% for each of DHFR12 and DHFR3. It was 300% for DB and 60% for AD in BFG-Y2H. DB\n",
    "dip_pos = sim.positive_interaction(dip,0.10,CV_autoactivity,pos_rate,cells_for_dox_induction) #CV of X-Y pair-dependent growth amplitudes (log-normal distribution) in BFG-PCA-selectable condition is 10%\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "prime-breeding",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Examine the coverage for diploids \n",
    "cov_pos = sim.dip_coverage(dip_pos,run_name)\n",
    "cov_pos_ppi = sim.dip_coverage([strain for strain in dip_pos if (strain[0][2]==1)],run_name)\n",
    "cov_neg = sim.dip_coverage(dip_neg,run_name)\n",
    "\n",
    "print(\"The ratio of diploid strains having at least the number of cells carried over to next step above threshold.\")\n",
    "print(\"\\nNon-selection:\")\n",
    "sim.print_LL(cov_neg)\n",
    "\n",
    "print(\"\\nSelection :\")\n",
    "sim.print_LL(cov_pos)\n",
    "\n",
    "print(\"\\nCoverage of strains with positives PPIs in selectable condition\")\n",
    "sim.print_LL(cov_pos_ppi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "focal-slovakia",
   "metadata": {},
   "outputs": [],
   "source": [
    "# You will be able to plot a heatmap of abundance at this stage."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "otherwise-healing",
   "metadata": {},
   "outputs": [],
   "source": [
    "###################################################\n",
    "# Plasmid prep\n",
    "\n",
    "\n",
    "plasmid_molecules = 1.4 * 10**7 #See below for estimation\n",
    "\n",
    "###################################################\n",
    "# Molecule estimation when performing 4 PCR reactions after plasmid prep;\n",
    "#\n",
    "# 100 pg / µL * 16  µL * 4 reactions = 6400 pg\n",
    "# Fraction of Y2H plasmids in yeast DNA miniprep product is 6% of the total DNA mass\n",
    "# Barcode fusion efficiency is 20%\n",
    "# 6400 pg * 0.06 * 0.2  = 152.6 pg\n",
    "# Y2H plasmid sizes are 10 kbp\n",
    "# 1 mole of double-stranded DNA molecule is 660 g bp^-1 (1ng = 0.152 fmol)\n",
    "# 153.6 * (1/1000 * 0.152) fmol  = 153.6 * 10**-3 * 0.152 * 10**-15 * 6.03* 10**23\n",
    "#                                = 153.6 * 6.03 * 0.152 * 10**5\n",
    "#                                = 1.4 * 10**7 molecules\n",
    "######################################################\n",
    "\n",
    "\n",
    "plasmid_neg = sim.adjust_num(dip_neg,0.5,plasmid_molecules)\n",
    "plasmid_pos = sim.adjust_num(dip_pos,0.5,plasmid_molecules)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "integral-request",
   "metadata": {},
   "outputs": [],
   "source": [
    "###############################################################\n",
    "# PCR re-amplification and Illumina sequencing\n",
    "\n",
    "reads_per_condition = 2* 10**6 # Reads for each condition (add Barcode Fusion type (BC1-BC1 or BC2-BC2) in each condition.)\n",
    "\n",
    "\n",
    "read_neg = sim.adjust_num(plasmid_neg,0.5,reads_per_condition)\n",
    "read_pos = sim.adjust_num(plasmid_pos,0.5,reads_per_condition)\n",
    "read_pos_PPI = sim.adjust_num([strain for strain in read_pos if (strain[0][2]==1)],0.5,reads_per_condition)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "increasing-appendix",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Examining coverage at final output\n",
    "cov_pos = sim.dip_coverage(read_pos,run_name)\n",
    "cov_pos_ppi = sim.dip_coverage([strain for strain in read_pos if (strain[0][2]==1)],run_name)\n",
    "cov_neg = sim.dip_coverage(read_neg,run_name)\n",
    "\n",
    "print(\"The ratio of diploid strains having at least the number of reads sequenced.\")\n",
    "print(\"\\nNon-selection:\")\n",
    "sim.print_LL(cov_neg)\n",
    "\n",
    "print(\"\\nSelection :\")\n",
    "sim.print_LL(cov_pos)\n",
    "\n",
    "print(\"\\nCoverage of strains with positives PPIs in selectable condition\")\n",
    "sim.print_LL(cov_pos_ppi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "abroad-shopper",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Some data formatting functions to enable data represenation\n",
    "\n",
    "# Haploid abundance histogram\n",
    "marginal_x_hap,marginal_y_hap = sim.marginal(read_neg)\n",
    "# Diploid abundance histogram\n",
    "abundance_pos                 = sim.abundance(read_neg)\n",
    "abundance_neg                 = sim.abundance(read_pos)\n",
    "\n",
    "# Function to format array data for heatmap representation in each condition \n",
    "hmap_read_neg = sim.format_array(read_neg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "referenced-habitat",
   "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.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}