notebooks/process_barcoderatios.py.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "77d9b03c-8057-4671-89c1-e5a883195895",
   "metadata": {},
   "source": [
    "# Process plate counts to get ratios of variants and determine pooling and MOI\n",
    "This notebook is designed to be run interactively, to determine the relative concentration of strains in equal volume and repooled samples\n",
    "\n",
    "## Setup\n",
    "Import Python modules:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "900fe170-1b4b-4646-b9ee-01d75a76f664",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pickle\n",
    "import sys\n",
    "\n",
    "import altair as alt\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import numpy\n",
    "\n",
    "import pandas as pd\n",
    "from os.path import join\n",
    "import os\n",
    "import ruamel.yaml as yaml\n",
    "\n",
    "_ = alt.data_transformers.disable_max_rows()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ecc47332-2e29-4605-bf8b-463edb63cac6",
   "metadata": {},
   "source": [
    "Get the variables passed by `snakemake`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "08570ac2-fd9d-42eb-8a19-8c62ecb27e97",
   "metadata": {},
   "outputs": [],
   "source": [
    "viral_library_csv = snakemake.input.viral_library\n",
    "neut_standard_set_csv = snakemake.input.neut_standard_set\n",
    "initialpool_metadata = snakemake.input.initialpool_metadata\n",
    "initial_pool_counts = snakemake.input.initial_pool_counts\n",
    "initial_pool_fates = snakemake.input.initial_pool_fates\n",
    "repooledlibraryfile = snakemake.input.repooledlibraryfile\n",
    "strainrepooling_volumes = snakemake.output.strainpooling\n",
    "equalvolume_strainbalance_plot = snakemake.output.equalvolume_plot\n",
    "repool_strainbalance_plot = snakemake.output.repool_plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "57502565-9d0d-4e85-a095-00b3b10d1af6",
   "metadata": {},
   "outputs": [],
   "source": [
    "barcode_runs_df = pd.read_csv(initialpool_metadata)\n",
    "barcode_runs_df.drop(columns=[\"fastq\"], inplace=True)\n",
    "barcode_runs_df[\"sample\"] = barcode_runs_df.apply(\n",
    "    lambda x: \"-\".join(x.astype(str)), axis=1\n",
    ")\n",
    "\n",
    "samples = barcode_runs_df[\"sample\"].unique().tolist()\n",
    "print(f\"There are {len(samples)} barcode runs.\")\n",
    "\n",
    "count_csvs = initial_pool_counts\n",
    "\n",
    "fate_csvs = initial_pool_fates\n",
    "\n",
    "samples_df = barcode_runs_df"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6bc6f5cd-de3d-4824-a421-7b460090d72e",
   "metadata": {},
   "source": [
    "## Statistics on barcode-parsing for each sample\n",
    "Make interactive chart of the \"fates\" of the sequencing reads parsed for each sample on the plate.\n",
    "\n",
    "If most sequencing reads are not \"valid barcodes\", this could potentially indicate some problem in the sequencing or barcode set you are parsing.\n",
    "\n",
    "Potential fates are:\n",
    " - *valid barcode*: barcode that matches a known virus or neutralization standard, we hope most reads are this.\n",
    " - *invalid barcode*: a barcode with proper flanking sequences, but does not match a known virus or neutralization standard. If you  have a lot of reads of this type, it is probably a good idea to look at the invalid barcode CSVs (in the `./results/barcode_invalid/` subdirectory created by the pipeline) to see what these invalid barcodes are.\n",
    " - *unparseable barcode*: could not parse a barcode from this read as there was not a sequence of the correct length with the appropriate flanking sequence.\n",
    " - *low quality barcode*: low-quality or `N` nucleotides in barcode, could indicate problem with sequencing.\n",
    " - *failed chastity filter*: reads that failed the Illumina chastity filter, if these are reported in the FASTQ (they may not be).\n",
    "\n",
    "Also, if the number of reads per sample is very uneven, that could indicate that you did not do a good job of balancing the different samples in the Illumina sequencing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9aa103c7-3240-4dac-8532-f7298a5fe009",
   "metadata": {},
   "outputs": [],
   "source": [
    "fates = (\n",
    "    pd.concat([pd.read_csv(f).assign(well=os.path.basename(f).strip('_fates.csv')) for f, s in zip(fate_csvs, samples)])\n",
    "    .merge(samples_df, validate=\"many_to_one\", on=\"well\")\n",
    "    .assign(\n",
    "        fate_counts=lambda x: x.groupby(\"fate\")[\"count\"].transform(\"sum\"),\n",
    "        sample_well=lambda x: x[\"sample\"] + \" (\" + x[\"well\"] + \")\",\n",
    "    )\n",
    "    .query(\"fate_counts > 0\")[  # only keep fates with at least one count\n",
    "        [\"fate\", \"count\", \"well\", \"sample_well\", \"dilution_factor\"]\n",
    "    ]\n",
    ")\n",
    "\n",
    "assert len(fates) == len(fates.drop_duplicates())\n",
    "\n",
    "\n",
    "sample_wells = list(\n",
    "    fates.sort_values([\"dilution_factor\"])[\"sample_well\"]\n",
    ")\n",
    "\n",
    "\n",
    "\n",
    "fates_chart = (\n",
    "    alt.Chart(fates)\n",
    "    .encode(\n",
    "        alt.X(\"count\", scale=alt.Scale(nice=False, padding=3)),\n",
    "        alt.Y(\n",
    "            \"sample_well\",\n",
    "            title=None,\n",
    "            sort=sample_wells,\n",
    "        ),\n",
    "        alt.Color(\"fate\", sort=sorted(fates[\"fate\"].unique(), reverse=True)),\n",
    "        alt.Order(\"fate\", sort=\"descending\"),\n",
    "        tooltip=fates.columns.tolist(),\n",
    "    )\n",
    "    .mark_bar(height={\"band\": 0.85})\n",
    "    .properties(\n",
    "        height=alt.Step(10),\n",
    "        width=200,\n",
    "        title=f\"Barcode parsing for initial titering plate\",\n",
    "    )\n",
    "    .configure_axis(grid=False)\n",
    ")\n",
    "\n",
    "fates_chart"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a810fb9a-189d-48e5-b866-6ef0f2e0dd5f",
   "metadata": {},
   "source": [
    "## Read barcode counts\n",
    "Read the counts per barcode:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ccea55e6-d280-4e8e-998a-1a86930c6529",
   "metadata": {},
   "outputs": [],
   "source": [
    "# get barcode counts\n",
    "counts = (\n",
    "    pd.concat([pd.read_csv(c).assign(well=os.path.basename(c).strip('_counts.csv')) for c, s in zip(count_csvs, samples)])\n",
    "    .merge(samples_df, validate=\"many_to_one\", on=\"well\")\n",
    "    .drop(columns=[\"replicate\"])\n",
    "    .assign(sample_well=lambda x: x[\"sample\"] + \" (\" + x[\"well\"] + \")\")\n",
    ")\n",
    "\n",
    "# classify barcodes as viral or neut standard\n",
    "barcode_class = pd.concat(\n",
    "    [\n",
    "        pd.read_csv(viral_library_csv)[[\"barcode\", \"strain\"]].assign(\n",
    "            neut_standard=False,\n",
    "        ),\n",
    "        pd.read_csv(neut_standard_set_csv)[[\"barcode\"]].assign(\n",
    "            neut_standard=True,\n",
    "            strain=pd.NA,\n",
    "        ),\n",
    "    ],\n",
    "    ignore_index=True,\n",
    ")\n",
    "\n",
    "# merge counts and classification of barcodes\n",
    "assert set(counts[\"barcode\"]) == set(barcode_class[\"barcode\"])\n",
    "counts = counts.merge(barcode_class, on=\"barcode\", validate=\"many_to_one\")\n",
    "assert set(sample_wells) == set(counts[\"sample_well\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e1326e88-a654-4a21-b9e5-85f5bf93847f",
   "metadata": {},
   "source": [
    "Apply any manually specified data drops:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "faedb6a2-ebcd-4d56-8ce9-dddfc2a0569d",
   "metadata": {},
   "source": [
    "## Average counts per barcode in each well"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7ddb2544-a1a3-460f-9d81-fa061b152071",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-12-05T19:26:10.449100Z",
     "iopub.status.busy": "2023-12-05T19:26:10.448471Z",
     "iopub.status.idle": "2023-12-05T19:26:10.468023Z",
     "shell.execute_reply": "2023-12-05T19:26:10.467117Z",
     "shell.execute_reply.started": "2023-12-05T19:26:10.449041Z"
    }
   },
   "source": [
    "Plot average counts per barcode.\n",
    "If a sample has inadequate barcode counts, it may not have good enough statistics for accurate analysis, and a QC-threshold is applied:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "68da99a3-e1e4-4b03-a7c2-79406e6d7170",
   "metadata": {},
   "outputs": [],
   "source": [
    "avg_barcode_counts = (\n",
    "    counts.groupby(\n",
    "        [\"well\", \"sample_well\"],\n",
    "        dropna=False,\n",
    "        as_index=False,\n",
    "    )\n",
    "    .aggregate(avg_count=pd.NamedAgg(\"count\", \"mean\"))\n",
    "    .assign(\n",
    "        fails_qc=lambda x: (\n",
    "            x[\"avg_count\"] < 500\n",
    "        ),\n",
    "    )\n",
    ")\n",
    "\n",
    "avg_barcode_counts_chart = (\n",
    "    alt.Chart(avg_barcode_counts)\n",
    "    .encode(\n",
    "        alt.X(\n",
    "            \"avg_count\",\n",
    "            title=\"average barcode counts per well\",\n",
    "            scale=alt.Scale(nice=False, padding=3),\n",
    "        ),\n",
    "        alt.Y(\"sample_well\", sort=sample_wells),\n",
    "        alt.Color(\n",
    "            \"fails_qc\",\n",
    "            title=f\"fails {'min barcode count threshold'=}\",\n",
    "            legend=alt.Legend(titleLimit=500),\n",
    "        ),\n",
    "        tooltip=[\n",
    "            alt.Tooltip(c, format=\".3g\") if avg_barcode_counts[c].dtype == float else c\n",
    "            for c in avg_barcode_counts.columns\n",
    "        ],\n",
    "    )\n",
    "    .mark_bar(height={\"band\": 0.85})\n",
    "    .properties(\n",
    "        height=alt.Step(10),\n",
    "        width=250,\n",
    "        title=f\"Average barcode counts per well for titering plate\",\n",
    "    )\n",
    "    .configure_axis(grid=False)\n",
    ")\n",
    "\n",
    "display(avg_barcode_counts_chart)\n",
    "\n",
    "# drop wells failing QC\n",
    "avg_barcode_counts_per_well_drops = list(avg_barcode_counts.query(\"fails_qc\")[\"well\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "589a755d-b46f-46ea-a84c-2dd784fa8672",
   "metadata": {},
   "source": [
    "## Fraction of counts from neutralization standard\n",
    "Determine the fraction of counts from the neutralization standard in each sample, and make sure this fraction passess the QC threshold."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "29a71183-4670-40c1-88c1-237bb536803e",
   "metadata": {},
   "outputs": [],
   "source": [
    "neut_standard_fracs = (\n",
    "    counts.assign(\n",
    "        neut_standard_count=lambda x: x[\"count\"] * x[\"neut_standard\"].astype(int)\n",
    "    )\n",
    "    .groupby(\n",
    "        [\"well\", \"sample_well\"],\n",
    "        dropna=False,\n",
    "        as_index=False,\n",
    "    )\n",
    "    .aggregate(\n",
    "        total_count=pd.NamedAgg(\"count\", \"sum\"),\n",
    "        neut_standard_count=pd.NamedAgg(\"neut_standard_count\", \"sum\"),\n",
    "    )\n",
    "    .assign(\n",
    "        neut_standard_frac=lambda x: x[\"neut_standard_count\"] / x[\"total_count\"],\n",
    "        fails_qc=lambda x: (\n",
    "            x[\"neut_standard_frac\"] < 0.001\n",
    "        ),\n",
    "    )\n",
    ")\n",
    "\n",
    "neut_standard_fracs_chart = (\n",
    "    alt.Chart(neut_standard_fracs)\n",
    "    .encode(\n",
    "        alt.X(\n",
    "            \"neut_standard_frac\",\n",
    "            title=\"frac counts from neutralization standard per well\",\n",
    "            scale=alt.Scale(nice=False, padding=3),\n",
    "        ),\n",
    "        alt.Y(\"sample_well\", sort=sample_wells),\n",
    "        alt.Color(\n",
    "            \"fails_qc\",\n",
    "            title=f\"fails {'min_neut_standard_frac_per_well'=}\",\n",
    "            legend=alt.Legend(titleLimit=500),\n",
    "        ),\n",
    "        tooltip=[\n",
    "            alt.Tooltip(c, format=\".3g\") if neut_standard_fracs[c].dtype == float else c\n",
    "            for c in neut_standard_fracs.columns\n",
    "        ],\n",
    "    )\n",
    "    .mark_bar(height={\"band\": 0.85})\n",
    "    .properties(\n",
    "        height=alt.Step(10),\n",
    "        width=250,\n",
    "        title=f\"Neutralization-standard fracs per well for titering plate\",\n",
    "    )\n",
    "    .configure_axis(grid=False)\n",
    "    .configure_legend(titleLimit=1000)\n",
    ")\n",
    "\n",
    "display(neut_standard_fracs_chart)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9b976631-4e94-4324-9ba7-46628ac09b8d",
   "metadata": {},
   "outputs": [],
   "source": [
    "counts_balancedbarcode = counts.groupby(['sample','strain','dilution_factor','serum','well','barcode']).sum()\n",
    "counts_balancedbarcode = counts_balancedbarcode.reset_index().drop(columns = ['neut_standard'])\n",
    "counts_balancedbarcode"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d940310f-9774-4508-b504-55d0eacf4a22",
   "metadata": {},
   "source": [
    "## Fraction of counts from each barcode\n",
    "Determine the fraction of counts from the each barcode for in each well of the plate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4ae846d6-71f6-4252-9a4a-008085d25a4c",
   "metadata": {},
   "outputs": [],
   "source": [
    "#First calculate the total counts per barcode in each well\n",
    "sumperwell = counts_balancedbarcode.groupby(['sample','sample_well','dilution_factor','serum','well']).sum().drop(columns=['strain','barcode','date'])\n",
    "sumperwell = sumperwell.reset_index()\n",
    "sumperwell"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "93c09a44-2aac-4bab-b9b5-bd3cb65f96ea",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Then calculate the sum of barcodes corresponding to barcoded variants per well and calculate the fraction of counts that each barcode represents of barcode counts per well\n",
    "sumperwell = sumperwell.rename(columns={'count':'counts_perwell'})\n",
    "counts_balancedbarcode = counts_balancedbarcode.merge(sumperwell, on=['sample','sample_well','dilution_factor','serum','well'])\n",
    "counts_balancedbarcode['fraction_barcode'] = counts_balancedbarcode['count'] /counts_balancedbarcode['counts_perwell'] \n",
    "counts_balancedbarcode"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2e017d45-d521-4a01-a96f-3e25b55572b7",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Now we select a couple wells where we think the vRNA counts should correspond to linearly to titers based on prior experiments amd calculate the sum of all barcodes for each strain and determine the fraction of counts that is each strain\n",
    "selected_well = counts_balancedbarcode.loc[counts_balancedbarcode['sample'].str.contains('-A6-')] #Take wells where fraction neut-standard is reasonable\n",
    "sum_barcodes_bystrain = selected_well.groupby(['strain','well'])['fraction_barcode'].sum().to_frame().rename(columns = {'fraction_barcode': 'fraction_strain_perwell'}).reset_index()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "77a6a161-6b83-4ca3-8247-84568483a371",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Merge calculation back with other data on barcoded strains\n",
    "mean_single_well = selected_well.merge(sum_barcodes_bystrain, on = ['strain','well'], how = 'left')\n",
    "\n",
    "# calcualte ratio to add\n",
    "num_strains = 36\n",
    "mean_single_well['ratio_to_add'] = (1/num_strains)/mean_single_well['fraction_strain_perwell']\n",
    "\n",
    "mean_single_well"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "98c5ddd6-c154-4d30-a046-78c26b128500",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#Plot the fraction that each strain is represented\n",
    "(mean_single_well\n",
    "                 [['strain', 'fraction_strain_perwell']]\n",
    "                 .drop_duplicates()\n",
    "                 .plot.barh(x=\"strain\", y=\"fraction_strain_perwell\", figsize = (6,10), log=False, xlim =(0,0.1))\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9f8ba2fb-14c5-42a5-a6c6-7899c649c430",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#Make a dataframe that just has the fraction each strain represented and the ratio is should be added to the pool\n",
    "initial_pool_ratios = mean_single_well[['strain', 'fraction_strain_perwell', 'ratio_to_add']].drop_duplicates()\n",
    "initial_pool_ratios.plot.barh(x=\"strain\", y=\"ratio_to_add\", figsize = (6,14))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "da8ce3de-9269-4790-882d-64bc536414c0",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Calculate volume to add to each well\n",
    "initial_pool_ratios['vol_to_add'] = initial_pool_ratios['ratio_to_add'].apply(lambda x: x * 200).to_list()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ea6a47bb-7df1-4b3f-b7d3-f7f6b52becdf",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Send volumes to a printable csv file\n",
    "initial_pool_ratios[['strain','vol_to_add']].to_csv(strainrepooling_volumes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b96ce4cf-0d5f-442e-81d9-96a4f7e3ea71",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Look at the representation of each barcode in this well\n",
    "mean_single_well['strain_barcode'] = mean_single_well['strain'] + \"_\" + mean_single_well['barcode'].str[0:16]\n",
    "mean_single_well.plot.barh(x=\"strain_barcode\", y=\"fraction_barcode\", figsize = (6,40), log=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "03c3b62f-ed57-4ce2-bba9-96b5257a449e",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Now we want to generate a plot that shows barcodes by strain\n",
    "assess_barcodebalancing = mean_single_well[['strain','barcode','strain_barcode','fraction_barcode','count','counts_perwell']]\n",
    "assess_barcodebalancing_bystrain = assess_barcodebalancing.groupby(['strain'])['count'].sum().to_frame().rename(columns = {'count': 'count_perstrain_perwell'}).reset_index()\n",
    "assess_barcodebalancing_bystrain\n",
    "assess_barcodebalancing = assess_barcodebalancing.merge(assess_barcodebalancing_bystrain, on=\"strain\")\n",
    "assess_barcodebalancing['fraction_of_strain_barcode'] = assess_barcodebalancing['count']/assess_barcodebalancing['count_perstrain_perwell']\n",
    "assess_barcodebalancing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "82676861-3668-48cd-b608-abdd47cff245",
   "metadata": {},
   "outputs": [],
   "source": [
    "barcode_balance_chart = (\n",
    "    alt.Chart(assess_barcodebalancing)\n",
    "    .encode(\n",
    "        alt.X(\"fraction_of_strain_barcode\", scale=alt.Scale(nice=False, padding=3)),\n",
    "        alt.Y(\n",
    "            \"strain\",\n",
    "            title=None,\n",
    "            sort=sample_wells,\n",
    "        ),\n",
    "        alt.Color(\"barcode\", sort=sorted(assess_barcodebalancing[\"strain\"].unique())).scale(range = ['steelblue', 'goldenrod','firebrick','rebeccapurple']),\n",
    "    )\n",
    "    .mark_bar(height={\"band\": 0.75})\n",
    "    .properties(\n",
    "        height=alt.Step(20),\n",
    "        width=250,\n",
    "    )\n",
    "    .configure_axis(grid=False,labelFontSize=15,titleFontSize=18,labelLimit=300)\n",
    ")\n",
    "\n",
    "barcode_balance_chart\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ee8adc57-2417-4461-a209-d92128cd7223",
   "metadata": {},
   "source": [
    "# Confirm that the ratio of the strains represented is similar after repooling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d6aff639-6aa1-4182-8456-bf5f1a3452f9",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Checking how well repooling worked based on balancing in plate 1\n",
    "counts_repooled = repooledlibraryfile\n",
    "counts_repooled = pd.read_csv(counts_repooled)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "028fa7cf-87a9-4364-a323-158e5b393a7e",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Merge with strain names\n",
    "counts_repooled_df = counts_repooled.merge(barcode_class, on=\"barcode\", validate=\"many_to_one\")\n",
    "counts_repooled_df = counts_repooled_df.loc[counts_repooled_df['neut_standard'] == False].drop(columns=['neut_standard'])\n",
    "counts_repooled_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a901ca66-274e-4508-9b14-d5371c93d262",
   "metadata": {},
   "outputs": [],
   "source": [
    "counts_repooled_df_bystrain = counts_repooled_df.groupby(['strain']).sum(numeric_only=True).reset_index()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b6e04150-c261-4903-a514-44d63dddd218",
   "metadata": {},
   "outputs": [],
   "source": [
    "counts_repooled_df_bystrain['counts_well'] = counts_repooled_df_bystrain['count'].sum()\n",
    "\n",
    "counts_repooled_df_bystrain['fraction'] = counts_repooled_df_bystrain['count']/counts_repooled_df_bystrain['counts_well']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f2b739df-08ba-46ef-bd27-b70b5a2d1544",
   "metadata": {},
   "outputs": [],
   "source": [
    "counts_repooled_df_bystrain.plot.barh(x=\"strain\", y=\"fraction\", figsize = (6,10), xlim=(0,0.1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "321864f2-470a-41d7-b0b2-f53a8a1233b5",
   "metadata": {},
   "outputs": [],
   "source": [
    "#We can look at this also as difference from the desired fraction, rather than raw fraction\n",
    "counts_repooled_df_bystrain['ratio_to_desired_fraction'] = counts_repooled_df_bystrain['fraction']/(1/36)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "37d16900-d469-44b2-bdfd-59c34503c50e",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Also calculate ratio to desired fraction in original pool\n",
    "initial_pool_ratios['ratio_to_desired_fraction'] = initial_pool_ratios['fraction_strain_perwell']/(1/36)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "46905a52-80cd-47aa-a737-135c5177afc6",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Plot difference in representation in original pool and in repooled sample\n",
    "counts_repooled_df_bystrain.plot.barh(x=\"strain\", y=\"ratio_to_desired_fraction\", xlim=(0,3), figsize = (4,10))\n",
    "initial_pool_ratios.plot.barh(x=\"strain\", y=\"ratio_to_desired_fraction\", xlim=(0,3),figsize = (4,10))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "201e741d-5d0f-4b6d-8c7b-647dc9d9041f",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Plot difference in representation in original pool and in repooled sample\n",
    "counts_repooled_df_bystrain.plot.barh(x=\"strain\", y=\"fraction\", xlim=(0,0.1), figsize = (4,10))\n",
    "plt.savefig(repool_strainbalance_plot,\n",
    "            dpi = 'figure',\n",
    "            bbox_inches = 'tight')\n",
    "\n",
    "initial_pool_ratios.plot.barh(x=\"strain\", y=\"fraction_strain_perwell\", xlim=(0,0.1),figsize = (4,10))\n",
    "plt.savefig(equalvolume_strainbalance_plot,\n",
    "            dpi = 'figure',\n",
    "            bbox_inches = 'tight')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "99fc8623-240a-4769-9a5e-32afcc5ef947",
   "metadata": {},
   "source": [
    "# Determine concentration of library that should be used for neutralization assay"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "41591442-db68-4cbf-8557-ed4f95dfda83",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Determine when we are in the linear range for vRNA relationship to sample dilutions\n",
    "MOItestsamples = neut_standard_fracs.loc[neut_standard_fracs['sample_well'].str.contains('230417')]\n",
    "MOItestsamples['TCID50'] = MOItestsamples['sample_well'].str.split('-').str[3].str.split(\"_\").str[1].str[:-1]\n",
    "MOItestsamples.plot.scatter(x='TCID50', y='neut_standard_frac')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "354199e4-080f-477e-a7b0-98a234ad57b2",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}