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",
    "import numpy as np\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": [
    "# Defining 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": [
    "# Defning number of haploid cells subjected to mating\n",
    "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",
    "cov = sim.hap_coverage(x_hap,y_hap,run_name)"
   ]
  },
  {
   "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": "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
}