template_notebook.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Template `Jupyter` notebook for computing stability scores from experimental data uploaded to TACC\n",
    "\n",
    "This notebook is a template for computing stability scores starting from raw FACS and deep-sequencing data uploaded by BIOFAB to TACC.\n",
    "\n",
    "There are five main sections:\n",
    "\n",
    "1. Define all input variables\n",
    "2. Compile all metadata into a single `experiments.csv` file for use as input for computing EC50 values\n",
    "3. Compute EC50 values from the experimental data\n",
    "4. Compute stabiltiy scores from EC50 values\n",
    "5. Generate summary plots\n",
    "\n",
    "As an example, this notebook analyzes data from a previous experiment with data located within the `sd2e-community/` directory on TACC:\n",
    "\n",
    "    sd2e-community/shared-q1-workshop/strcklnd/protein_design_data/Plan_18338\n",
    "\n",
    "The `sd2e-community/` directory is accessible to anyone. For instance, it is located at the followig path for me:\n",
    "    \n",
    "    /work/05402/haddox/jupyter/\n",
    "    \n",
    "Hugh Haddox, October-8-2018"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import `Python` modules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import sys\n",
    "sys.path.append('scripts/')\n",
    "import prot_stab_utils\n",
    "import re\n",
    "import glob\n",
    "\n",
    "import xml.etree.ElementTree as ET\n",
    "sys.path.append(\"/home/jupyter/tacc-work/jupyter_packages/lib/python2.7/site-packages\")\n",
    "from FlowCytometryTools import *\n",
    "\n",
    "import numpy as np\n",
    "import math\n",
    "import scipy.stats\n",
    "import pandas\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "import seaborn as sns\n",
    "sns.set(font_scale=2)\n",
    "sns.set_style(\"ticks\")\n",
    "from IPython.display import Image"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1) Define all input variables\n",
    "\n",
    "The below cell contains ***all*** input variables for the notebook. Just change these variables and then run the below cells in an automated fashion.\n",
    "\n",
    "### Input variables for compiling experimental metadata\n",
    "\n",
    "Typically, the BIOFAB uploads all data in a single folder, which in turn contains the following files/sub-directories:\n",
    "\n",
    "* `manifest.csv`: a file that lists all samples tested in the experiment, along with metadata for each sample, such as unique IDs for locating files with associated FACS and deep-sequencing data.\n",
    "* `ngs_data/`: a directory with all deep-sequencing data.\n",
    "* `facs_data/`: a directory with all FACS data.\n",
    "\n",
    "Below are input variables that are relevant to compiling these data:\n",
    "\n",
    "* `summary_file`: a path to a CSV file with metadata for each sample\n",
    "* `library_name`: the name of the library in the `strain` column of the `summary_file`. There are usually multiple values, including the library of interest, as well as controls (e.g., `AMA1-best`).\n",
    "\n",
    "* `map_job_id_to_sort_round`: a dictionary that is used to match values in the `sort_job` column of `manifest.csv` (dictionary key) to the corresponding round of sorting (dictionary value). For instance, the original experiment from Rocklin et al., 2017, Science, involves three rounds of sorting: the first round with the naive library and the libraries treated with the two lowest concentrations of protease, the second round with library treated with the two intermediate concentrations of protease, and the third round with the library treated with the two highest concentrations of protease.\n",
    "\n",
    "* `sort_rounds_with_identical_control`: a string of comma-delimited integers giving sorting rounds where the same sample was used as the protease-naive control for both proteases. These samples provide a baseline for how many cells are expected to pass the sorting threshold in the *absence* of protease. Typically, this is true for the first round of sorting, in which case `1` would be an appropriate entry.\n",
    "\n",
    "* `fastq_dir`: a path to a folder with the deep-sequencing data\n",
    "\n",
    "* `facs_dir`: a path to a folder with the FACS data\n",
    "\n",
    "* `resultsdir`: the path to a directory where all results from this script will be stored. This directory will be created if it does not already exist.\n",
    "\n",
    "* `experimental_metadata_output_file`: the path of the output file that will be created. This file is similar to the `experiments.csv` file from Rocklin et al, 2017, Science."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define input files and directories\n",
    "data_dir = '/work/05402/haddox/jupyter/sd2e-community/shared-q1-workshop/strcklnd/protein_design_data/Plan_18338'\n",
    "summary_file = os.path.join(data_dir, 'manifest.csv')\n",
    "fastq_dir = os.path.join(data_dir, 'ngs_data/')\n",
    "facs_dir = os.path.join(data_dir, 'facs_data/')\n",
    "\n",
    "# Define variables that are relevant for parsing the\n",
    "# input summary file\n",
    "library_name = 'protein design downselect 3 with ladder 1'\n",
    "map_job_id_to_sort_round = {\n",
    "    'Job_72019' : 1, 'Job_72263' : 2, 'Job_72471' : 3\n",
    "}\n",
    "sort_rounds_with_identical_control = [1]\n",
    "\n",
    "# Define the path to the output summary file\n",
    "resultsdir = \"results/template_notebook\"\n",
    "experimental_metadata_output_file = os.path.join(resultsdir, 'experimental_metadata_from_script.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Input variables for computing EC50 values and stability scores from experimental data\n",
    "\n",
    "* `designed_sequences_file` : the path to a CSV file giving the name and protein sequence for each input design. See [here](data/Rocklin_2017_Science/designed_protein_sequences.csv) for an example. In this file, each row specifies a design and each column specifies information about that design. This file must have the following comma-delimited columns:\n",
    "    * `name` : a unique name for the design\n",
    "    * `protein_sequence` : the protein sequence of the design\n",
    "\n",
    "* `pare_path` : a path to the program [`PEAR`](https://sco.h-its.org/exelixis/web/software/pear/)\n",
    "\n",
    "* `five_prime_flanking_seq`: a DNA sequence that flanks the coding sequence of interest on the 5' end of the sequencing read (string). The coding sequence should begin immediately after the final nucleotide of this flanking sequence. This flanking sequence and the one given by `three_prime_flanking_seq` will be used to extract the DNA coding sequence from each sequencing read and then translate it into a protein dsequence. Note: the default sequence used in Rocklin et al., 2017, Science was `CATATG`.\n",
    "\n",
    "* `three_prime_flanking_seq`: a DNA sequence that flanks the coding sequence of interest on the 3' end of the sequencing read (string). The coding sequence should begin immediately before the first nucleotide of this flanking sequence. Note: this DNA sequence should be in the same 5'-to-3' orientation as `five_prime_flanking_seq`. Note: the default sequence used in Rocklin et al., 2017, Science was `CTCGAG`\n",
    "\n",
    "* `conc_factor`: the fold-change in protease concentration between selection steps, i.e., when `selection_strength` is incrimented by a value of one. A value of 3 would indicate that the protease concentration is increased by 3 fold between selection steps. Note: leave this value blank for samples that have not been treated with any protease.\n",
    "\n",
    "* `proteases`: a list of the proteases used in the experiment that are lower-case versions of the names in the file specified by the `summary_file` input variable.\n",
    "\n",
    "* `sorting_threshold`: the fluorescence threshold (log10 space) used to sort cells as \"positive\" or \"negative\" in the FACS step of the experiment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "designed_sequences_file = 'data/template_notebook/designed_sequences.csv'\n",
    "pear_path = '/home/05402/haddox/software/pear/bin/pear'\n",
    "five_prime_flanking_seq = 'CATATG'\n",
    "three_prime_flanking_seq = 'CTCGAG'\n",
    "conc_factor = '3'\n",
    "proteases = ['trypsin', 'chymotrypsin']\n",
    "sorting_threshold = np.log10(2000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2) Compile all metadata into a single `experiments.csv` file for use as input for computing EC50 values\n",
    "\n",
    "The goal of this section is to compile all metadata on TACC uploaded by the BIOFAB for a given experiment. The result of this section is an `experiments.csv` file, similar to the one from Rocklin et al., that is used as input for computing EC50 values. For more details on this file, see the documentation for the script `compute_ec50_values_from_deep_sequencing_data.py` in the `README.md` file at the root of this repository."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### First, read in metadata from the input `manifest.csv` file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read in data\n",
    "summary_df = pandas.read_csv(summary_file)\n",
    "\n",
    "# The manifest usually has an entry for a protease-naive sample.\n",
    "# Sometimes the BIOFAB leaves all entries for this sample blank,\n",
    "# except for the `aq_item_id` column. In case this is true for this\n",
    "# manifest, I replace any nan entries with naive in the protease\n",
    "# column, also filling in the strain column. I make sure that there\n",
    "# is only one naive sample at the end. If the naive sample already\n",
    "# had an entry in the manifest, this shouldn't change anything\n",
    "summary_df['protease'] = summary_df['protease'].fillna('naive')\n",
    "summary_df['strain'] = summary_df['strain'].fillna(library_name)\n",
    "assert sum(summary_df['protease'] == 'naive') == 1\n",
    "\n",
    "# For each sample, specify a unique file prefix that will be used to find FASTQ files for that sample\n",
    "summary_df['aq_item_id'] = summary_df['aq_item_id'].fillna(0).astype(int)\n",
    "summary_df['fastq_id'] = summary_df['aq_item_id'].apply(\n",
    "    lambda x: os.path.join(str(x), 'Files/')\n",
    ")\n",
    "\n",
    "# For each sample, specify a unique file prefix that will be used to find FACS files for that sample\n",
    "summary_df.rename(columns={'sort_job':'job_id'}, inplace=True)\n",
    "summary_df['facs_file_prefix'] = summary_df.apply(\n",
    "    lambda row: os.path.join(facs_dir, str(row['job_id']), str(row['facs_filename_stub'])), axis=1\n",
    ")\n",
    "\n",
    "# Make other changes to the dataframe\n",
    "summary_df['protease'] = summary_df['protease'].apply(lambda x: x.lower())\n",
    "summary_df.set_index(['protease', 'concentration'], inplace=True)\n",
    "summary_df.sort_index(inplace=True)\n",
    "\n",
    "# Remove entries that do not correspond to the library of interest\n",
    "summary_df = summary_df[\n",
    "    summary_df['strain'] == library_name\n",
    "]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Analyze the FACS data, recording relevant information for each sample such as the fraction of cells passing the FITC cutoff for both the sample and untreated controls, as well as the total number of cells collected for each sample"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_facs_file(facs_file_prefix):\n",
    "    \"\"\"Find all files with a given prefix\"\"\"\n",
    "    \n",
    "    # Find all files with given prefix, make sure there\n",
    "    # is only one matching file, and return its path\n",
    "    facs_files = glob.glob(facs_file_prefix + '*.fcs')\n",
    "    if len(facs_files) == 1:\n",
    "        return facs_files[0]\n",
    "    elif len(facs_files) == 0:\n",
    "        return None\n",
    "    \n",
    "def s_log(sample, channel_names):\n",
    "    \"\"\"Transform FACS data to log10 space\"\"\"\n",
    "    \n",
    "    # Copy dataframe and get relevant FACS data\n",
    "    new_sample = sample.copy()\n",
    "    new_data = new_sample.data\n",
    "\n",
    "    # Transform the FACS data to log10 space\n",
    "    for channel_name in channel_names:\n",
    "        new_data[channel_name] = np.log10(new_data[channel_name])\n",
    "    \n",
    "    # Replace values of negative infinity with negative one, remove\n",
    "    # entries with `nan` values, and return a dataframe with the\n",
    "    # final values\n",
    "    new_data.replace(to_replace = -np.inf, value = -1, inplace = True)\n",
    "    new_data = new_data.dropna()\n",
    "    new_sample.data = new_data\n",
    "    \n",
    "    return new_sample\n",
    "\n",
    "def compute_fraction_collected(fcs_file, sorting_threshold=3.0):\n",
    "    \"\"\"Compute the fraction of cells that passed the FITC gate during sorting\"\"\"\n",
    "    \n",
    "    # Return `None` if the file doesn't exist\n",
    "    if not isinstance(fcs_file, str):\n",
    "        return None\n",
    "    elif not os.path.isfile(fcs_file):\n",
    "        return None\n",
    "    \n",
    "    # Get FACS data from a given file and for a given chanel\n",
    "    else:\n",
    "        sample = FCMeasurement(ID='Test Sample', datafile=fcs_file)\n",
    "        sample = s_log(sample, [facs_channel_name])\n",
    "        df = sample.data\n",
    "\n",
    "        # Count the number of total events and events passing the gate\n",
    "        total_number_of_recorded_events = float(len(df.index.values))\n",
    "        number_of_events_passing_fitc_gate = sum(df[facs_channel_name] > sorting_threshold)\n",
    "        fraction_collected = number_of_events_passing_fitc_gate / total_number_of_recorded_events\n",
    "        #return (pandas.Series([total_number_of_recorded_events, number_of_events_passing_fitc_gate, fraction_collected]))\n",
    "        return fraction_collected"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, compute the fraction of cells that were collected for each sample"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Search for the FACS file associated with each sample\n",
    "summary_df['facs_file_name'] = summary_df['facs_file_prefix'].apply(\n",
    "    lambda x: find_facs_file(x)\n",
    ")\n",
    "\n",
    "# Specify the FACS channel name to investigate\n",
    "facs_channel_name = u'FITC-A'\n",
    "\n",
    "# Compute the faction of cells that passed the FITC gate for each sample\n",
    "summary_df['fraction_collected'] = summary_df.apply(\n",
    "    lambda row: compute_fraction_collected(\n",
    "        row['facs_file_name'],\n",
    "        sorting_threshold\n",
    "    ), axis=1\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the distribution of FITC values for each sample with this type of data, as well as the sorting threshold used to select cells"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make a dataframe with only samples that have FACS data\n",
    "# and then initiate a figure with subplots\n",
    "data = summary_df[~summary_df['facs_file_name'].isnull()].copy()\n",
    "nplots = len(data)\n",
    "ncols = 3\n",
    "nrows = int(math.ceil(nplots/float(ncols)))\n",
    "xticks = list(range(0,7))\n",
    "fig = plt.figure(figsize=[15,20])\n",
    "\n",
    "# Plot the FACS data, showing the sorting threshold in each\n",
    "# case\n",
    "for (plot_n, (i, row)) in enumerate(data.iterrows(), 1):\n",
    "\n",
    "    # Get metadata and read in FACS data\n",
    "    (protease_type, concentration) = i\n",
    "    fcs_file = row['facs_file_name']\n",
    "    selection_round = map_job_id_to_sort_round[row['job_id']]\n",
    "    sample = FCMeasurement(ID='Test Sample', datafile=fcs_file)\n",
    "    sample = s_log(sample, [facs_channel_name])\n",
    "    df = sample.data\n",
    "    \n",
    "    # Plot the FACS data\n",
    "    ax = fig.add_subplot(nrows, ncols, plot_n)\n",
    "    sns.distplot(df[facs_channel_name], ax=ax)\n",
    "    sns.despine()\n",
    "    ax.set_title('{0}_{1}_{2}'.format(protease_type, concentration, selection_round))\n",
    "    ax.set_xlabel('log10 FITC')\n",
    "    ax.set_ylabel('density of cells')\n",
    "    ax.plot([sorting_threshold, sorting_threshold], [0, 1.0], c='k', ls='--')\n",
    "axs = fig.get_axes()\n",
    "plt.setp(axs, xticks=xticks)\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute `parent_expression` values. For each sample, this value gives the fraction of cells passing the FITC cutoff in a protease-naive untreated control sample, setting a baseline for how many cells would be expected to pass in the absence of protease. After computing these values, I will remove the naive samples used to compute these values as they are not used in further analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make a dictionary that gives `parent_expression` values for each sorting round\n",
    "print(\"The following rounds have identical naive controls: {0}\".format(\n",
    "    ', '.join(map(str, sort_rounds_with_identical_control))\n",
    "))\n",
    "naive_sample_expression = {'trypsin':{}, 'chymotrypsin':{}}\n",
    "list_of_indices_to_drop = []\n",
    "\n",
    "for (i, row) in summary_df.iterrows():\n",
    "    (protease, concentration) = i\n",
    "    if math.isnan(concentration):\n",
    "        continue\n",
    "    sort_round = map_job_id_to_sort_round[row['job_id']]\n",
    "    if concentration == 0.0:\n",
    "        assert (sort_round not in naive_sample_expression[protease].keys()), \"Found multiple controls for same round\"\n",
    "        \n",
    "        # If a single control is used for both proteases, record an identical value\n",
    "        # for both proteases. Otherwise, record protease-specific values\n",
    "        if sort_round in sort_rounds_with_identical_control:\n",
    "            for protease_i in naive_sample_expression.keys():\n",
    "                naive_sample_expression[protease_i][sort_round] = row['fraction_collected']\n",
    "        else:\n",
    "            naive_sample_expression[protease][sort_round] = row['fraction_collected']\n",
    "        list_of_indices_to_drop.append(i)\n",
    "        \n",
    "# Remove protease naive samples that will not be used in downstream steps\n",
    "summary_df.drop(list_of_indices_to_drop, axis='index', inplace=True)\n",
    "summary_df.reset_index(inplace=True)\n",
    "\n",
    "# Add the `parent_expression` column\n",
    "def record_parent_expression(protease, job_id):\n",
    "    if job_id not in map_job_id_to_sort_round:\n",
    "        return None\n",
    "    else:\n",
    "        sort_round = map_job_id_to_sort_round[job_id]\n",
    "        if sort_round in naive_sample_expression[protease].keys():\n",
    "            return naive_sample_expression[protease][sort_round]\n",
    "        else:\n",
    "            return None\n",
    "        \n",
    "summary_df['parent_expression'] = summary_df.apply(\n",
    "    lambda row: record_parent_expression(row['protease'], row['job_id']), axis=1\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, I will determine the total number of cells collected for each sample. These data are encoded in separate files than the ones analyzed above. In these files, the data are organized by specimen and tube numbers. These numbers are also present in the above dataframe in the column called `facs_filename_stub`. I will parse these values and then use them to match the new FACS data with the appropriate samples in the dataframe."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def parse_specimen_and_tube(spec_and_tube_string):\n",
    "    \"\"\"Parse the specimen and tube numbers from an input string with an expected pattern\"\"\"\n",
    "    \n",
    "    # Define the expected pattern\n",
    "    spec_and_tube_pattern = re.compile(r'(?P<spec>Specimen_\\d+)_(?P<tube>Tube_\\d+)')\n",
    "    \n",
    "    # Return `None` if the input variable isn't a string or does not\n",
    "    # contain the expected pattern. Otherwise, return the parsed values\n",
    "    if not isinstance(spec_and_tube_string, str):\n",
    "        return(pandas.Series([None, None]))\n",
    "    else:\n",
    "        match = re.search(spec_and_tube_pattern, spec_and_tube_string)\n",
    "        if match:\n",
    "            return(pandas.Series([match.group('spec'), match.group('tube')]))\n",
    "        else:\n",
    "            return(pandas.Series([None, None]))\n",
    "\n",
    "summary_df[['specimen', 'tube']] = summary_df.apply(\n",
    "    lambda row: parse_specimen_and_tube(row['facs_filename_stub']), axis=1\n",
    ")\n",
    "\n",
    "summary_df.set_index(['job_id', 'specimen', 'tube'], inplace=True)\n",
    "summary_df.sort_index(inplace=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Read in data from the XML files on the total number of cells collected. I will do this for all XML files the BIOFAB has uploaded, even ones that aren't relevant to this experiment. Then, I will identify the relevant ones using the `job_id`, `specimen`, and `tube` columns."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Compile data in each XML file\n",
    "xml_data_dict = {\n",
    "    key : []\n",
    "    for key in ['job_id', 'specimen', 'tube', 'cells_collected', 'intended_number_of_cells_collected'] #\n",
    "}\n",
    "cells_collected_pattern = re.compile(r'\\w+ : (?P<cells_collected>\\d+) / (?P<intended_number_of_cells_collected>\\d+)')\n",
    "xml_files = glob.glob(os.path.join(facs_dir, '*/*/*.xml'))\n",
    "for filename in xml_files:\n",
    "    tree = ET.parse(filename)\n",
    "    root = tree.getroot()\n",
    "\n",
    "    # Get experiment metadata\n",
    "    xml_data_dict['job_id'].append( root.findall(\".//*[@name='Experiment']\")[0].get('value') )\n",
    "    xml_data_dict['specimen'].append( root.findall(\".//*[@name='Specimen']\")[0].get('value') )\n",
    "    xml_data_dict['tube'].append( root.findall(\".//*[@name='Tube']\")[0].get('value') )\n",
    "    xml_data_dict['xmlfile'] = filename\n",
    "    \n",
    "    # Get experiment counts\n",
    "    cells_collected_data = root[4].findall(\".//*[@col='1']\")[0].text\n",
    "    match = re.match(cells_collected_pattern, cells_collected_data)\n",
    "    if match:\n",
    "        xml_data_dict['cells_collected'].append(\n",
    "            int(match.group('cells_collected'))\n",
    "        )\n",
    "        xml_data_dict['intended_number_of_cells_collected'].append(\n",
    "            int(match.group('intended_number_of_cells_collected'))\n",
    "        )\n",
    "    else:\n",
    "        xml_data_dict['cells_collected'].append(None)\n",
    "        xml_data_dict['intended_number_of_cells_collected'].append(None)\n",
    "          \n",
    "# Convert data to a dataframe\n",
    "xml_data_df = pandas.DataFrame.from_dict(xml_data_dict)\n",
    "xml_data_df.set_index(['job_id', 'specimen', 'tube'], inplace=True)\n",
    "\n",
    "# Check for entries with duplicate indices and remove the one with the\n",
    "# fewest counts if multiple exist\n",
    "if sum(xml_data_df.index.duplicated()) > 0:\n",
    "    print(\"Found duplicate indices. Will keep the duplicate with the highest counts, removing the rest\")\n",
    "    xml_data_df.sort_values(by='cells_collected', ascending=False, inplace=True)\n",
    "    xml_data_df = xml_data_df[~xml_data_df.index.duplicated()]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Merge the FACS data from the two input sources"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Merge the two sources of FACS data\n",
    "summary_df = summary_df.merge(\n",
    "    xml_data_df[['cells_collected', 'intended_number_of_cells_collected']],\n",
    "    left_index=True, right_index=True, how='left'\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Tweak the names of certain entries and manually add missing columns"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Change the entry for the protease-naive unsorted sample to have the `protease` column equal `trypsin` and the `concentration` column to equal `0.0`, which is consistent with the naming scheme used in Rocklin et al., 2017, Science."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add a row for a naive sample for each protease\n",
    "trypsin_naive_df = summary_df[\n",
    "    summary_df['protease'] == 'naive'\n",
    "].replace('naive', 'trypsin')\n",
    "\n",
    "chymotrypsin_naive_df = summary_df[\n",
    "    summary_df['protease'] == 'naive'\n",
    "].replace('naive', 'chymotrypsin')\n",
    "summary_df = pandas.concat([summary_df, trypsin_naive_df, chymotrypsin_naive_df])\n",
    "\n",
    "# Convert the concentration of the naive samples to `0.0`, drop the row with `naive`\n",
    "# in the `protease` column, then sort by protease and concentration\n",
    "summary_df['concentration'].replace(np.nan, 0.0, inplace=True)\n",
    "summary_df.reset_index(inplace=True)\n",
    "summary_df.set_index(['protease', 'concentration'], inplace=True)\n",
    "summary_df.drop(('naive', 0.0), inplace=True)\n",
    "summary_df.sort_index(inplace=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Manually add missing columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "summary_df['experiment_id'] = library_name\n",
    "summary_df['selection_strength'] = 2 * [i for i in range(0,7)]\n",
    "summary_df['conc_factor'] = 2 * ['', '3', '3', '3', '3', '3', '3']\n",
    "summary_df['parent'] = 2 * ['', '0', '0', '2', '2', '4', '4']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Rename columns, write the final dataframe to an output file, and show the final dataframe."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Rename columns\n",
    "summary_df.reset_index(inplace=True)\n",
    "summary_df['protease_type']=summary_df['protease']\n",
    "\n",
    "column_order = [\n",
    "    'experiment_id', 'protease_type', 'concentration', 'selection_strength', 'parent', 'conc_factor', 'fastq_id',\n",
    "    'parent_expression', 'fraction_collected', 'cells_collected'\n",
    "]\n",
    "\n",
    "# Write the dataframe to an output file\n",
    "if not os.path.isdir(resultsdir):\n",
    "    os.makedirs(resultsdir)\n",
    "print(\"Writing metadata to the experimental summary file: {0}\".format(experimental_metadata_output_file))\n",
    "summary_df[column_order].to_csv(experimental_metadata_output_file, index=False)\n",
    "\n",
    "# Show the final dataframe\n",
    "summary_df[column_order]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check that the `fraction_collected` column, computed above using the raw FACS data and the reported selection threshold, is close to the fraction recorded by the person who conducted the experiment, which is given in the `frac_positive` column."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "summary_df['facs_fractions_match'] = summary_df.apply(\n",
    "    lambda row: np.isclose(\n",
    "        row['frac_positive'], row['fraction_collected'],\n",
    "        atol=0.01, equal_nan=True\n",
    "    ), axis=1\n",
    ")\n",
    "if sum(~summary_df['facs_fractions_match']) > 0:\n",
    "    print(\"FACS values do not match for the following entries\")\n",
    "    display(summary_df[\n",
    "        ~summary_df['facs_fractions_match']\n",
    "    ][[\n",
    "        'protease_type', 'concentration', 'frac_positive',\n",
    "        'fraction_collected', 'facs_fractions_match'\n",
    "    ]])\n",
    "    raise ValueError()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## 3) Compute EC50 values from the experimental data\n",
    "\n",
    "In the below cell, I assemble and print a command-line argument for running the script `compute_ec50_values_from_deep_sequencing_data.py`.\n",
    "\n",
    "To actually run this script requires activating the `Conda` environment specified by the file called `environment_compute_ec50_values.yml` and then executing the arguement, which I will do manually."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Write the command to carry everything out\n",
    "cmd = ' '.join([\n",
    "    'python',\n",
    "    'scripts/compute_ec50_values_from_deep_sequencing_data.py',\n",
    "    '--designed_sequences_file {0}'.format(designed_sequences_file),\n",
    "    '--experimental_summary_file {0}'.format(experimental_metadata_output_file),\n",
    "    '--fastq_dir {0}'.format(fastq_dir),\n",
    "    '--pear_path {0}'.format(pear_path),\n",
    "    '--five_prime_flanking_seq {0}'.format(five_prime_flanking_seq),\n",
    "    '--three_prime_flanking_seq {0}'.format(three_prime_flanking_seq),\n",
    "    '--output_dir {0}'.format(resultsdir)\n",
    "])\n",
    "print(\"To complete this step, manually activate the environment encoded in environment_compute_ec50_values.yml:\\n\")\n",
    "print('source activate 2018_prot_stab_compute_ec50_values')\n",
    "print(\"\\nThen run the command:\\n\")\n",
    "print(cmd)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4) Compute stability scores from EC50 values\n",
    "\n",
    "In the below cell, I assemble and print a command-line argument for running the script `compute_stability_scores_from_EC50_values.py`.\n",
    "\n",
    "To actually run this script requires activating the `Conda` environment specified by the file called `environment_compute_stability_scores.yml` and then executing the below arguements, which I will do manually."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Put together the command\n",
    "print(\"To complete this step, manually activate the environment encoded in environment_compute_stability_scores.yml\\n\")\n",
    "print('source activate 2018_prot_stab_compute_stability_scores')\n",
    "print(\"\\nThen run the commands:\\n\")\n",
    "for protease in proteases:\n",
    "\n",
    "    # Define the names of the input file with EC50 values and the output file\n",
    "    # that will have the computed stability scores\n",
    "    ec50_values_file = os.path.join(\n",
    "        resultsdir,\n",
    "        'ec50_values/{0}.fulloutput'.format(protease)\n",
    "    )\n",
    "    output_file = os.path.join(\n",
    "        resultsdir,\n",
    "        'stability_scores/{0}_stability_scores.txt'.format(protease)\n",
    "    )\n",
    "\n",
    "    # Define and print the command\n",
    "    cmd = ' '.join([\n",
    "        'python',\n",
    "        'scripts/compute_stability_scores_from_EC50_values.py',\n",
    "        designed_sequences_file,\n",
    "        protease,\n",
    "        ec50_values_file,\n",
    "        conc_factor,\n",
    "        output_file\n",
    "    ])\n",
    "\n",
    "    print(cmd)\n",
    "    print('\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, merge stability scores between proteases, computing the minimum stability scores for the two proteases, and creating a file with a single dataframe with data from each protease, with suffixes of `_t` or `_c` indicating that the data is derived from trypsin or chymotrypsin, respectively, and a new column called `stabilityscore` that gives the minimum stability score between the two proteases."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stability_scores_trypsin_file = os.path.join(\n",
    "    resultsdir,\n",
    "    'stability_scores/trypsin_stability_scores.txt'\n",
    ")\n",
    "stability_scores_chymotrypsin_file = os.path.join(\n",
    "    resultsdir,\n",
    "    'stability_scores/chymotrypsin_stability_scores.txt'\n",
    ")\n",
    "merged_stability_scores_outfile = os.path.join(\n",
    "    resultsdir,\n",
    "    'stability_scores/stability_scores.txt'\n",
    ")\n",
    "print(\"Writing merged stability scores to the file: {0}\".format(merged_stability_scores_outfile))\n",
    "prot_stab_utils.merge_stability_scores_between_proteases(\n",
    "    stability_scores_trypsin_file,\n",
    "    stability_scores_chymotrypsin_file,\n",
    "    merged_stability_scores_outfile\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5) Make summary plots\n",
    "\n",
    "In the below cell, I define input variables and then assemble and print a command-line argument for running the script `create_summary_plots.py`.\n",
    "\n",
    "To actually run this script requires activating the `Conda` environment specified by the file called `environment_compute_ec50_values.yml` and then executing the arguement, which I will do manually."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Define input variables\n",
    "data_dir = resultsdir\n",
    "output_dir = os.path.join(resultsdir, 'summary_plots/')\n",
    "\n",
    "# Write the command to carry everything out\n",
    "cmd = ' '.join([\n",
    "    'python',\n",
    "    'scripts/create_summary_plots.py',\n",
    "    '--data_dir {0}'.format(data_dir),\n",
    "    '--output_dir {0}'.format(output_dir)\n",
    "])\n",
    "\n",
    "print(\"To complete this step, manually activate the environment encoded in environment_compute_ec50_values.yml:\\n\")\n",
    "print('source activate 2018_prot_stab_compute_ec50_values')\n",
    "print(\"\\nThen run the command:\\n\")\n",
    "print(cmd)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Show the number of assembled paired-end reads per sample, as well as the number of reads that could not be assembled or were discarded."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Image(os.path.join(output_dir, \"deep_sequencing_depth_and_quality.png\"), width=700)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}