<jupyter_start><jupyter_text>Kaz Lab Python Plotting Workshop***2018.05.25Updated 2018.07.24 ObjectivesBy the end of this workshop, you will be able to:- Create, launch, and run commands in an interactive Python ([IPython](https://ipython.readthedocs.io/en/stable/index.html)) notebook- Load single timepoint or time series data from a plate reader using `pandas`- Convert plate reader data into the "tidy" format, and join your raw data to metadata in a plate map- Make exploratory plots which compare different variables and subsets of your data Un-objectivesThis workshop will **not** teach you Python in general; our goal is to learn just enough Python that you can make useful plots. Likewise we are going to use several giant and complicated libraries—but only very small parts of them.I will focus on quickly generating plots for exploratory data analysis—I will not cover the numerous options for generating publication- (or even presentation-) quality plots.  Getting started  The Notebook and YouLet's first get acquainted with the notebook. This notebook is part of a project called [Jupyter](http://jupyter.org), which, in general, lets you write and run code "interactively"---that is, you write a little snippet of code, then run it and see the results. You can keep repeating this, and your variables live on across multiple snippets. Behind the scenes, your notebook is sending these commands to a program—called a "kernel"—which executes the commands, keeps track of your variables, and sends you the results. (In this case, the kernel is [IPython](https://ipython.readthedocs.io/en/stable/index.html), but [you can use many other programming languages as well](https://github.com/jupyter/jupyter/wiki/Jupyter-kernels)!) **This means your variables stick around as long as the kernel stays alive.** If you close this browser window and return to it later, your variables will still be there. If you restart your computer, the notebook server, or the kernel, you will have to re-run the cells in your notebook to get your variables back.Try running the code in the cell below: click on the box, then press `Shift + Enter`<jupyter_code>1 + 1<jupyter_output><empty_output><jupyter_text>The output appears below the input cell. This highlights another feature of the notebook—**it is a record of what commands you ran to analyze your data**. In this way, it's a lot like a lab notebook for analysis. This is very important for a few reasons:1. **It makes your analysis reproducible**; given your notebook and your data, someone else could repeat the exact same analysis to verify your results.2. **It makes it analysis verifiable**; since there is a record of the exact commands you ran, someone else can read your code to see if you have made a mistake, an invalid assumption, etc. If you are copying and pasting data from Excel into Prism, there is no way for someone else to check that every range was copied correctly. If you make a mistake, it can be silently propagated through your thinking and analysis, and your conclusions could be completely wrong!  The very basics of Python:- Define a variable `x` and assign it a value like this:        ```python    x = 10    ```        Note that a single equals sign (`=`) is **assignment**—it _sets_ `x` equal to `10`. To test if     `x` _is already equal_ to `10`, use two equals signs (`==`):        ```ipython     In[1]: x == 10    Out[1]: True    ```        - Text surrounded by quotes (`'` or `"`) is called a string:    ```python    y = 'Hello'    ```        - You can call functions by writing the name of the function, followed by parenthesis, with any parameters (also called arguments) inside the parentheses—like this:    ```python    print("Hello world")    ```    - Anything that follows a `` sign is a comment:    ```python    print("Hello world")  stuff back here doesn't do anything    ``` A few more tips about the Notebook:- To get help with an object or function, type it's name followed by `?`. Try it in the box below:<jupyter_code>int?<jupyter_output><empty_output><jupyter_text>- IPython can auto-complete many commands if you type a few characters and press `tab`. **Exercise:** Run the cell below, then see if you can figure out how to capitalize the text in `s` so it reads `Hello world` (Hint: you can call functions like this `s.upper()`):<jupyter_code>s = 'hello world'
s.upper()
# finish this command so it prints 'Hello world'
s.cap<jupyter_output><empty_output><jupyter_text>- Use the white buttons near the top to add, cut/paste, and move cells. Click the `+` button to add a new cell- If you start running a command that's taking too long, use the `Kernel` menu to `Interrupt`. If things get really messed up, you can `Restart` the kernel. Once you've done that, you may want to use the `Cell` menu to run all the cells above your current position, in order to revive all your variables.- You can export your notebook to several different formats (like HTML or PDF), if you want to share your analysis with someone who doesn't have Jupyter/IPython.  Tools of the TradeEnough of that, let's get started. First, we will **`import`** a bunch of libraries so we can work with them.<jupyter_code>import matplotlib.pyplot as plt
import pandas as pd
import seaborn as sns
import numpy as np
import datetime
import re

import plates

%matplotlib inline
%config InlineBackend.figure_format = 'retina'
sns.set_context("notebook", font_scale=1.2)
pd.set_option("display.max_rows",10)<jupyter_output><empty_output><jupyter_text>I copy and paste almost this exact statement at the top of each of my notebooks; when I need other things, I'll add stuff, but this imports most of the important libraries that I use regularly. It's important having a big-picture conceptual understanding of what each library does, since they build on each other—sometimes a function that doesn't exist in a "higher-level" library can be found in a "lower-level" library.- [NumPy](https://docs.scipy.org/doc/numpy/) (`np`): provides an N-dimensional array (`ndarray`)—much like a matrix in MATLAB. Your data is a matrix of numbers, and NumPy provides that matrix- [Pandas](http://pandas.pydata.org/pandas-docs/stable/) (`pd`): provides objects Series, DataFrame, and Panel that wrap NumPy arrays and label their contents. For instance, rather than having to remember that column `7` of your matrix represents the strain, you can label that column as such. This may seem small, but it makes a huge difference. - [Matplotlib](https://matplotlib.org/api/pyplot_summary.html) (`plt`): low-level plotting library, originally based on the plotting functions in MATLAB. - [Seaborn](https://seaborn.pydata.org/api.html) (`sns`): high-level plotting library for making complex visualizations—most importantly "grids" of plots (a.k.a. facets, lattice plots, trellis plots, or "small multiples"), and statistical plots.- `plates` is a library I wrote for describing plate maps and joining them to your data. I think it makes the analysis more straightforward, but you could do all of this without it  Loading and "Wrangling" dataIn this workshop, we'll focus on data from a microplate reader, like the Tecan. Other data will be slightly different, but the overall process will be the same. Furthermore, we'll "wrangle" our data into a common format that will let us approach many data sets in the same manner.The basic approach:- Import raw data using Pandas- Convert raw data into "tidy" format- Join raw data to metadata (plate map)We'll work through several "case studies" based on experiments done in the lab. These were chosen to emphasize different types of data and techniques for analyzing/plotting with them—not because of the significance of the experiments or the data.  Case 1: OxyraseGoal: The goal of this experiment (CG016e) was to determine if a difference in growth rate (manifested by a difference in OD600 after ~16 hours of incubation) could be measured in non-shaking cultures of _Pseudomonas_ growing in closed Eppendorf tubes. Method: I inoculated exponential-phase _Pseudomonas_ into aliquots of LB media supplemented with various combinations of nitrate, oxyrase, and tungstate. For each media condition, I did 3 biological replicates. The next day, I vortexed each cultures and divided that culture into 3 wells of a 96-well plate, then measured the OD with the Tecan. Data: single, endpoint measurement of a 96-well plate  Case 1: Loading & WranglingLet's start by importing the data from the Excel file produced by the Tecan<jupyter_code>data = pd.read_excel('./data/2018-03-30_CG016e_doctored.xlsx', # path to the Excel file
                     sheet_name=0, # index of the sheet to read; you could also pass a name like 'Sheet1'
                     header=26, # which row (zero-indexed) to start reading from; skips all the extra metadata
                                # junk that Tecan inserts
                     skipfooter=4, # how many rows to skip at the bottom of the file (Tecan puts some stuff 
                                   # there too)
                     index_col=0)  # which column to use to "index" the data
data<jupyter_output><empty_output><jupyter_text>Neat, it looks like a plate. Let's discuss the data structures we're using and introduce some terminology:- This dataset is now in a form called a [Pandas **DataFrame**](https://pandas.pydata.org/pandas-docs/stable/generated/pandas.DataFrame.html), which is a 2D table- Notice that the dataset is labeled—each row has a label, and each column has a label. The list of labels for the rows (`A`, `B`, `C`, etc.) is called the **`index`**. The labels for the columns is just called **`columns`**.- You can pull out a single column from the DataFrame like this, which produces a 1D table called a [Pandas Series](https://pandas.pydata.org/pandas-docs/stable/generated/pandas.Series.html):<jupyter_code>data[1]<jupyter_output><empty_output><jupyter_text>Let's talk about how the data is formatted. Each row and column in the DataFrame corresponds to a row and column in our physical microtiter plate; initially this seems intuitive, but it will be difficult to work with the data in this form for a few reasons:- **There's no metadata**: we can't tell from this dataset what cells correspond to what conditions, which cells are replicates of one another, and so on.- Even if I told you that `A1:A12` corresponded to the `LB alone` condition, with `A1:A3`, `A4:A6`, and `A7:A9` each three technical replicates of three biological replicates, it would be tricky to pull those nine values and group them appropriately for plotting or analysis. **The structure of our table does not correspond to the structure of our data**, even if it mimics the physical layout of our experiment. - It would be even more difficult if I told you that rows `B`, `D`, `F`, and `H` had oxyrase, and I wanted you to compare them to the rows which did not; comparing these rows would require a different set of strategies than comparing the replicates in the same row. **Different comparisons require very different code**.To address these problems, we're going to convert our data into a ["Tidy" format](http://vita.had.co.nz/papers/tidy-data.pdf), which will make it easier to process, visualize, and analyze.<jupyter_code># collect data into "tidy" format, with one row per observation, one column for each variable (well, time, OD600)
data.columns = data.columns.rename('column')
data.index = data.index.rename('row')
data = data.reset_index()
data = data.melt(id_vars=['row'],var_name='column',value_name='OD600')
data['well'] = data['row'] + data['column'].map(str)
data['row'] = data['row'].map(plates.letters2row)
data['column'] = data['column'] - 1
data<jupyter_output><empty_output><jupyter_text>(Note: For the purposes of this tutorial, I've told Pandas to [only show a few rows of each dataset](https://pandas.pydata.org/pandas-docs/stable/options.htmlfrequently-used-options) to keep things clean. Here are some functions to show more and to show all of a data set; if the output gets long enough, the notebook will collapse cell, and you can scroll separately inside it.)<jupyter_code>def display_more(df):
    with pd.option_context("display.max_rows",40,"display.max_columns", 40):
        display(df)

def display_all(df):
    with pd.option_context("display.max_rows",None,"display.max_columns", None):
        display(df)<jupyter_output><empty_output><jupyter_text>At first, this doesn't seem like an improvement. We'll really see the gains when we attach metadata from our platemap to this raw data. Here's the full platemap:| _     | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | Media | Nitrate | Tungstate | Oxyrase ||-------|---|---|---|---|---|---|---|---|---|----|----|----|-------|---------|-----------|---------|| **A** | 1 | 1 | 1 | 2 | 2 | 2 | 3 | 3 | 3 | s  | s  | s  | LB    | -       | -         | -       || **B** | 1 | 1 | 1 | 2 | 2 | 2 | 3 | 3 | 3 | s  | s  | s  | OLB   | -       | -         | +       || **C** | 1 | 1 | 1 | 2 | 2 | 2 | 3 | 3 | 3 | s  | s  | s  | LBN   | +       | -         | -       || **D** | 1 | 1 | 1 | 2 | 2 | 2 | 3 | 3 | 3 | s  | s  | s  | OLBN  | +       | -         | +       || **E** | 1 | 1 | 1 | 2 | 2 | 2 | 3 | 3 | 3 | s  | s  | s  | LBT   | -       | +         | -       || **F** | 1 | 1 | 1 | 2 | 2 | 2 | 3 | 3 | 3 | s  | s  | s  | OLBT  | -       | +         | +       || **G** | 1 | 1 | 1 | 2 | 2 | 2 | 3 | 3 | 3 | s  | s  | s  | LBNT  | +       | +         | -       || **H** | 1 | 1 | 1 | 2 | 2 | 2 | 3 | 3 | 3 | s  | s  | s  | OLBNT | +       | +         | +       |Numbers within the table (1, 2, 3) refer to biological replicates, and (s) refers to the sterile wells. Each row is a different media condition, with various combinations of nitrate, tungstate, and oxyrase. In the next function, we'll build a table which has the same structure as our raw data---one row for every well---but we'll add columns corresponding to our other "metadata" variables:<jupyter_code># build platemap with relevant metadata
platemap = plates.prog2spec({
  "A:A": {"condition":"LB",    "oxyrase": 0, "nitrate": 0,  "tungstate": 0},
  "B:B": {"condition":"OLB",   "oxyrase": 1, "nitrate": 0,  "tungstate": 0},
  "C:C": {"condition":"LBN",   "oxyrase": 0, "nitrate": 10, "tungstate": 0},
  "D:D": {"condition":"OLBN",  "oxyrase": 1, "nitrate": 10, "tungstate": 0},
  "E:E": {"condition":"LBT",   "oxyrase": 0, "nitrate": 0,  "tungstate": 10},
  "F:F": {"condition":"OLBT",  "oxyrase": 1, "nitrate": 0,  "tungstate": 10},
  "G:G": {"condition":"LBNT",  "oxyrase": 0, "nitrate": 10, "tungstate": 10},
  "H:H": {"condition":"OLBNT", "oxyrase": 1, "nitrate": 10, "tungstate": 10},
  "A1:H3": {"sample":"1"},
  "A4:H6": {"sample":"2"},      
  "A7:H9": {"sample":"3"}, 
  "A10:H12": {"sample":"sterile", "sterility":'sterile' }, 
}, include_row_column=False)
platemap<jupyter_output><empty_output><jupyter_text>All this confusing-looking statement does is turn a list of statements about our plate into a DataFrame which organizes those statements into a table---one row per well, one column for each type of metadata.  Finally, we'll merge our platemap and our raw data into a single table. The call to `pd.merge` tells Pandas to do an [inner join](https://en.wikipedia.org/wiki/Join_(SQL)Inner_join) (give rows that show up in both the `data` and `platemap` DataFrames), and to use the column `well` from `data`, and the `index` from the `platemap` to match rows between the two DataFrames<jupyter_code># join platemap to observations, 
data = pd.merge(data,platemap, how='inner', left_on='well', right_index=True)
data = data.fillna({'sterility': 'non-sterile'})

data<jupyter_output><empty_output><jupyter_text>Now we have our complete dataset in a "Tidy" format. **"Tidy" data** has the following form:- Each **variable** is in a column- Each **observation** is in a row. "A dataset is a collection of **values**, usually either numbers (if quantitative) or strings (if qualitative). Values are organised in two ways. Every value belongs to a variable and an observation. A **variable** contains all values that measure the same underlying attribute (like height, temperature, duration) across units. An **observation** contains all values measured on the same unit (like a person, or a day, or a race) across attributes" ([Wickam 2014](http://vita.had.co.nz/papers/tidy-data.pdf)).![](http://r4ds.had.co.nz/images/tidy-1.png)It's tricky to precisely define "observations" and "variables" in general, but it's usually easy to tell for a given dataset. In this example, our "observations" correspond to single wells, so each row contains all the data about a given well. Our variables are `OD600`, `well`, `oxyrase`, `nitrate`, `tungstate`, `sample`, and `sterility`. Notice that each is the name of a column in our table. Somewhat confusingly, we also have variables called `row` and `column`---these correspond to the row/column of each well within the original, physical microplate (Don't get confused between the rows/columns of wells in the physical plate and the rows/columns in our `DataFrame`!). We also defined a convenience column `condition`, which summarizes the information in the `oxyrase`, `nitrate`, and `tungstate` columns.   Case 1: Plotting For single timepoint data, I always like to start by looking at the data for the whole plate as a heatmap to see if any patterns jump out:<jupyter_code>g = sns.heatmap(data=data.pivot('row','column','OD600'))<jupyter_output><empty_output><jupyter_text>> (The `data.pivot` statement brings our data back into a form that looks a lot like the one we started with. It may seem silly to go through all this trouble to convert data into a tidy format—only to convert it straight back! But this demonstrates an important point—there are lots of functions built in to Pandas (and other libraries) to take tidy data and convert it into other formats.)<jupyter_code>data.pivot('row','column','OD600')<jupyter_output><empty_output><jupyter_text>It looks like the columns on the right have a much lower OD than the ones on the left—this makes sense, because  they were the sterile wells. Let's plot the mean of the sterile wells against the non-sterile wells:<jupyter_code>g = sns.catplot(data=data,x="sterility",y="OD600",kind="bar")<jupyter_output><empty_output><jupyter_text>[`catplot`](https://seaborn.pydata.org/generated/seaborn.catplot.html) is an incredibly useful function provided by Seaborn, which will be your Swiss-Army knife for [working with categorical data](https://seaborn.pydata.org/tutorial/categorical.html) (hence the name `catplot`, for `cat`egorical plot). Every call to `catplot` looks like this: - first pass it your dataset            g = sns.catplot(data=data,x="sterile",y="OD600",kind="bar")                        ^^^^^^^^^        - then tell it how to map your variables onto elements of the plot:        g = sns.catplot(data=data,x="sterile",y="OD600",kind="bar")                                  ^^^^^^^^^^^ ^^^^^^^^^- finally, tell it what `kind` of plot to draw        g = sns.catplot(data=data,x="sterile",y="OD600",kind="bar")                                                        ^^^^^^^^^^You can draw a surprising variety of [`kind`s](https://seaborn.pydata.org/api.htmlcategorical-plots) of plots, just with this one function. (note: `catplot` used to be called `factorplot`, so if you see references to `factorplot`, you can use it exactly the same way as `catplot`)Our bar plot collapses the entire distribution down to a mean and confidence interval. If we wanted to see the full range of values in the `'sterile'` and `'non-sterile'` groups, we could use a `box` plot:<jupyter_code>g = sns.catplot(data=data,x="sterility",y="OD600",kind="box")<jupyter_output><empty_output><jupyter_text>**Exercise** Try it yourself: adapt this code to make a `swarm`plot and a `violin`plot: What are the advantages of these plots over the `bar` plot or the `box` plots? ***Let's make a different comparison: are the oxyrase wells different from the non-oxyrase wells? **Exercise:** Generalize the code from above to make a `bar` plot with `oxyrase` concentration on the x-axis and `OD600` on the y-axis: Looks like they're not different in general, but what about if we separate the sterile vs. non-sterile groups:<jupyter_code>g = sns.catplot(data=data,x="sterility",y="OD600",hue="oxyrase",kind="bar")<jupyter_output><empty_output><jupyter_text>Notice we added the `hue` parameter in order to separate the color of bars based on the oxyrase concentration. Let's look at a `swarm` plot to see the distribution:<jupyter_code>g = sns.catplot(data=data,x="sterility",y="OD600",hue="oxyrase",kind="swarm")<jupyter_output><empty_output><jupyter_text>There are two populations within the `'non-sterile'` group. What could explain the difference between those two groups? Probably some combination of the `nitrate` and `tungstate` in the media. Let's break the plot into two plots, based on the nitrate concentration:<jupyter_code>g = sns.catplot(data=data,x="sterility",y="OD600",hue="oxyrase",col='nitrate',kind="swarm")<jupyter_output><empty_output><jupyter_text>This is where `catplot` really starts to shine—look how easily we were able to go from one plot to several small plots (or "facets"). We can go further too---let's make two rows of plots, one for each tungstate concentration:<jupyter_code>g = sns.catplot(data=data,x="sterility",y="OD600",hue="oxyrase",col='nitrate',row='tungstate',kind="swarm")<jupyter_output><empty_output><jupyter_text>What's more, you can easily generalize this code to make comparisons between different variables!**Exercise:** Switch it up, and adapt this code to plot `oxyrase` in rows, `sterility` in columns, `nitrate` on the x-axis, and `tungstate` concentration on different colors: Here are a bunch more example plots I made, using different combinations of variables:<jupyter_code>g = sns.catplot(data=data,hue="sample",x="condition",y="OD600",kind="swarm")
g.set_xticklabels(rotation=90)
g = sns.catplot(data=data,hue="sterility",x="condition",y="OD600",kind="bar")
g.set_xticklabels(rotation=90)
g = sns.catplot(data=data,hue="sterility",x="condition",y="OD600",kind="boxen",linewidth=0.5)
g.set_xticklabels(rotation=90)
g = sns.catplot(data=data,hue="sterility",x="condition",y="OD600",kind="violin",linewidth=0.5)
g.set_xticklabels(rotation=90)
g = sns.catplot(data=data,hue="sterility",x="condition",y="OD600",kind="box",linewidth=0.5)
g.set_xticklabels(rotation=90)
g = sns.catplot(data=data,hue="oxyrase",x="tungstate",y="OD600",col="nitrate", kind="box",linewidth=0.5)
g.set_xticklabels(rotation=90)<jupyter_output><empty_output><jupyter_text>To summarize what we learned about `catplot`: - use any combination of these these parameters to make plots with different comparisons:    - `x`: plot this variable on the x-axis; should be a categorical variable    - `y`: plot this variable on the y-axis; can be continuous or categorical    - `hue`: color elements of the plot according to this variable; categorical    - `col`: make columns of subplots according to this variable; categorical    - `row`: make rows of subplots according to this variable; categoricalChoose the `kind` of plot with `kind`. A [list of different categorical plots is here](https://seaborn.pydata.org/api.htmlcategorical-plots). Just remove `plot` from the end of the name (i.e. for `boxplot`, use `kind='box'`; for `boxenplot`, use `kind='boxen'`, etc.) ***One more point before we move on: We can also normalize the data to a certain condition. In this case, let's subtract the `sterile` replicate from the other replicates.<jupyter_code>import plates.data
data_normalized_sterile = plates.data.calc_norm(data,value='OD600',
                                    on='sterility',columns=['condition'], 
                                    how=lambda x: x - x.loc['sterile'].mean())

g = sns.catplot(data=data_normalized_sterile,hue="sample",x="condition",y="OD600",kind="swarm",
                   order=["LB","OLB","LBN","OLBN","LBT","OLBT","LBNT","OLBNT"])
g.set_xticklabels(rotation=90)<jupyter_output><empty_output><jupyter_text>Case 2: Amoeba co-incubationGoal: The goal of this experiment was to explore how co-incubation of two strains of _P. aeruginosa_---PA103 and PA14 FliF were affected by growth in the presence of _Acanthamoeba castelanii_ (*Ac*)Data: this experiment was performed over multiple 24-hour Tecan runs, then with several measurements taken irregularly so the data is in an Excel spreadsheet that contains several tables---one table for each set of reads.   Case 2: Loading & WranglingLet's start by building the plate map this time. The general layout of the plate is as follows:|     _      |   1    |   3    |   2    |   4    |   5    |   6    |     7     |     8     |     9     |     10    |     11    |     12    | *Ac* inoculum ||------------|--------|--------|--------|--------|--------|--------|-----------|-----------|-----------|-----------|-----------|-----------|---------------|| **Strain** | PA103  | PA103  | PA103  | PA103  | PA103  | PA103  | PA14 FliF | PA14 FliF | PA14 FliF | PA14 FliF | PA14 FliF | PA14 FliF |               || **A**      | $10^3$ | $10^3$ | $10^3$ | $10^4$ | $10^4$ | $10^4$ | $10^4$    | $10^4$    | $10^4$    | $10^3$    | $10^3$    | $10^3$    | $0$           || **B**      | $100$  | $100$  | $100$  | $10$   | $10$   | $10$   | $100$     | $100$     | $100$     | $10$      | $10$      | $10$      | $0$           || **C**      | $10^3$ | $10^3$ | $10^3$ | $10^4$ | $10^4$ | $10^4$ | $10^4$    | $10^4$    | $10^4$    | $10^3$    | $10^3$    | $10^3$    | $100$         || **D**      | $100$  | $100$  | $100$  | $10$   | $10$   | $10$   | $100$     | $100$     | $100$     | $10$      | $10$      | $10$      | $100$         || **E**      | $10^3$ | $10^3$ | $10^3$ | $10^4$ | $10^4$ | $10^4$ | $10^4$    | $10^4$    | $10^4$    | $10^3$    | $10^3$    | $10^3$    | $10^3$        || **F**      | $100$  | $100$  | $100$  | $10$   | $10$   | $10$   | $100$     | $100$     | $100$     | $10$      | $10$      | $10$      | $10^3$        || **G**      | $10^3$ | $10^3$ | $10^3$ | $10^4$ | $10^4$ | $10^4$ | $10^4$    | $10^4$    | $10^4$    | $10^3$    | $10^3$    | $10^3$    | $10^4$        || **H**      | $100$  | $100$  | $100$  | $10$   | $10$   | $10$   | $100$     | $100$     | $100$     | $10$      | $10$      | $10$      | $10^4$        |Numbers in the middle of the table are inoculum of _P. aeruginosa_. This platemap is a little more complicated than the first case. Check out the documentation of `plates.prog2spec` to see some of the possibilities for handling this plate:<jupyter_code>plates.prog2spec?
# build a platemap
platemap = plates.prog2spec({
    "1:6":  { "strain": "PA103" },
    "7:12": { "strain": "PA14 FliF" },
    "A:B":  { "Ac_inoculum": 0 },
    "C:D":  { "Ac_inoculum": 100 },
    "E:F":  { "Ac_inoculum": 10**3 },
    "G:H":  { "Ac_inoculum": 10**4 },
    "A:B,C:D,E:F,G:H": { "Pa_inoculum": [ 
        [ 10**3, 10**3, 10**3, 10**4, 10**4, 10**4, 10**4, 10**4, 10**4, 10**3, 10**3, 10**3 ],
        [ 100,   100,   100,   10,    10,    10,    100,   100,   100,   10,    10,    10,   ]
    ] }
}, include_row_column=True)

platemap<jupyter_output><empty_output><jupyter_text>Now we'll separately load each table from the Excel file, verifying that it contains the appropriate number of columns, and that it starts and ends at the right row.<jupyter_code>df1 = pd.read_excel('./data/5.5.18_growth curve_25C_PYG_121 hrs_Ac + PA103, PA14 FliF_MASTER.xlsx',
                    index_col=0,
                    header=2,
                    nrows=98,
                    usecols=plates.letters2row('EC'))
df1<jupyter_output><empty_output><jupyter_text>(I did this for each separate table, but for conciseness I am not showing it here)<jupyter_code>df2 = pd.read_excel('./data/5.5.18_growth curve_25C_PYG_121 hrs_Ac + PA103, PA14 FliF_MASTER.xlsx',
                    index_col=0,
                    header=104,
                    nrows=98,
                    usecols=plates.letters2row('EC'))
df3 = pd.read_excel('./data/5.5.18_growth curve_25C_PYG_121 hrs_Ac + PA103, PA14 FliF_MASTER.xlsx',
                    index_col=0,
                    header=206,
                    nrows=98,
                    usecols=plates.letters2row('CT'))
df4 = pd.read_excel('./data/5.5.18_growth curve_25C_PYG_121 hrs_Ac + PA103, PA14 FliF_MASTER.xlsx',
                    index_col=0,
                    header=309,
                    nrows=98,
                    usecols=3)
df5 = pd.read_excel('./data/5.5.18_growth curve_25C_PYG_121 hrs_Ac + PA103, PA14 FliF_MASTER.xlsx',
                    index_col=0,
                    header=412,
                    nrows=98,
                    usecols=3)
df6 = pd.read_excel('./data/5.5.18_growth curve_25C_PYG_121 hrs_Ac + PA103, PA14 FliF_MASTER.xlsx',
                    index_col=0,
                    header=514,
                    nrows=98,
                    usecols=4)<jupyter_output><empty_output><jupyter_text>This code looks complicated, but it's very similar to what we did for the first example: for each table, we'll convert into a tidy format, then merge with the platemap. There are a few more steps here, including extracting the temperature data and converting the times from seconds to hours<jupyter_code>offset = 0

# make a list of processed `DataFrame`s that will be combined together at the end
dfs = []

# make a list of temperatures which will be combined together at the end
temps = []

# for each of the consecutive 24-hour reads
for df in [df1, df2, df3]:
    
    # Extract temperature vs. time as separate variable, then drop (delete)
    # it from the main data frame
    temp = df.loc['Temp. [°C]',:]
    temp.index = (temp.index / 3600) + offset
    temps.append(temp)
    df = df.drop('Temp. [°C]')

    # likewise, extract the hours row that Carrie calculated;
    # we'll re-calculate it to check
    time_h = df.loc['Time [hrs]',:]
    df = df.drop('Time [hrs]')

    # convert columns to hours, and add offset
    df.columns = (df.columns / 3600) + offset
    
    # check that Carrie's calculation of the hours is correct
    assert np.allclose(df.columns.values,time_h.values)

    # the offset for the next table will be the last time point for this table
    offset = df.columns[-1]

    # collect df into "tidy" format, with one row per observation, 
    # one column for each variable (well, time, OD600)
    df.columns = df.columns.rename('time')
    df.index = df.index.rename('well')
    df = df.reset_index()
    df = df.melt(id_vars=['well'],value_name='OD600')
    
    # join platemap to the raw data
    df = pd.merge(df,platemap, how='inner', left_on='well',right_index=True)
    
    dfs.append(df)
    
# for each of the 
for df, offset in zip([df4, df5, df6], [88.5, 112, 121]):
    
    # Extract temperature vs. time as separate variable, then drop (delete)
    # it from the main data frame
    temp = df.loc['Temp. [°C]',:]
    temp.index = (temp.index / 3600) + offset
    temps.append(temp)
    df = df.drop('Temp. [°C]')
    
    # likewise, extract the hours row
    time_h = df.loc['Time [hrs]',:]
    df = df.drop('Time [hrs]')

    # convert columns to hours, and add offset
    df.columns = (df.columns / 3600) + offset

    # check that Carrie's calculation of the hours is correct
    assert np.allclose(df.columns.values,time_h.values)
    
    # collect df into "tidy" format, with one row per observation, 
    # one column for each variable (well, time, OD600)
    df.columns = df.columns.rename('time')
    df.index = df.index.rename('well')
    df = df.reset_index()
    df = df.melt(id_vars=['well'],value_name='OD600')
    
    # join platemap to the raw data
    df = pd.merge(df,platemap, how='inner', left_on='well',right_index=True)
    
    dfs.append(df)

temps = pd.concat(temps)
data = pd.concat(dfs)
    
display(data)
display(temps)<jupyter_output><empty_output><jupyter_text>We can compute new columns based on existing columns; let's calculate the multiplicity of infection (MOI):, as well as the log of the inocula:<jupyter_code>data['MOI'] = data['Pa_inoculum'] / data['Ac_inoculum']
data['log_MOI'] = np.log10(data['MOI'])
data['log_Pa_inoculum'] = np.log10(data['Pa_inoculum'])
data['log_Ac_inoculum'] = np.log10(data['Ac_inoculum']).apply(lambda x: 0 if x == float('-inf') else x)
data<jupyter_output><empty_output><jupyter_text>Case 2: Plotting The data in this example is different from the first case in a crucial way: the time axis is [continuous](https://en.wikipedia.org/wiki/Continuous_function), rather than [categorical (discrete)](https://en.wikipedia.org/wiki/Categorical_variable). Rather than using `catplot`, we'll meet a new function called [`relplot`](https://seaborn.pydata.org/generated/seaborn.relplot.htmlseaborn.relplot).Let's start by looking at a single condition: PA103, no *AC*, with an inoculum of 1000<jupyter_code>data.query("strain == 'PA103' & Ac_inoculum == 0 & Pa_inoculum == 1000")<jupyter_output><empty_output><jupyter_text>[`relplot`](https://seaborn.pydata.org/generated/seaborn.relplot.htmlseaborn.relplot) is the continuous cousin of `catplot`. Like `catplot`, we pass `relplot` our `data`, and then map figure aesthetics (such as `x`, `y`, and `hue`) to columns in our `data`. For example, let's map the time course of a single amoeba infection<jupyter_code>sns.relplot(data=data.query("strain == 'PA103' & Ac_inoculum == 0 & Pa_inoculum == 1000"),
             x='time',y='OD600',
             # stretch the plot a bit so we can see earlier timepoints more easily.
             aspect=1.5,
             # by default, Seaborn draws confidence intervals by Bootstrap sampling
             # <https://en.wikipedia.org/wiki/Bootstrapping_(statistics)>. This gets 
             # really slow with `relplot`, so we'll ask it to draw confidence intervals
             # using the standard deviation instead. To go back to bootstrapping, just remove
             # this line. It'll take much longer to plot!
             ci='sd',
             kind='line'
            )<jupyter_output><empty_output><jupyter_text>Turns out we can set `hue` on `relplot`, just like we did for `factorplot`; let's look at different starting inocula of PA103:<jupyter_code>sns.relplot(data=data.query("strain == 'PA103' & Ac_inoculum == 0"),
            x='time',y='OD600',hue='Pa_inoculum',
            aspect=1.5,
            palette='viridis',
            ci='sd',kind='line')<jupyter_output><empty_output><jupyter_text>This is probably not what you expected, but it makes sense when you realize that in `relplot`, everything is continuous. Under the hood, Seaborn is mapping values of your `hue` variable (`Pa_inoculum`) to a *continuous* linear scale. This turns out to be a very useful feature when the variable mapped to `hue` has lots of values---rather than showing each level of the hue as a separate entry in the legend, it can show a sampling of values.However, in this case it would be more helpful to scale the hues used to represent values of `Pa_inoculum` to a log scale, so that there's just as clear a difference between 1000 and 100 as between $10^4$ and 1000.<jupyter_code>from matplotlib.colors import LogNorm
sns.relplot(data=data.query("strain == 'PA103' & Ac_inoculum == 0"),
            x='time',y='OD600',hue='Pa_inoculum', 
            # normalize the colors to the log of the inoculum
            hue_norm=LogNorm(),
            aspect = 1.5,
            palette='viridis', 
            # hide the legend, we'll draw it ourself later
            legend=False,
            ci='sd',kind='line')

# override the legend with a nicer display of the inocula
plt.gca().legend(title='Pa inoculum (cfu)',labels=['$10$','$100$','$10^3$','$10^4$']);<jupyter_output><empty_output><jupyter_text>There are other ["normalization functions" for colormaps](https://matplotlib.org/users/colormapnorms.html) provided by [matplotlib](https://matplotlib.org).  `relplot` has a few other tricks up its sleeve. You can draw spaghetti plots like this:<jupyter_code>from matplotlib.colors import LogNorm
sns.relplot(data=data.query("strain == 'PA103' & Ac_inoculum == 0 & time < 20"),
            x='time',y='OD600',hue='Pa_inoculum',     
            hue_norm=LogNorm(),
            
            # plot each well as a distinct sampling unit
            units='well',
            # don't sum over wells
            estimator = None,

            aspect = 1.5,
            palette='viridis', 
            legend=False,
            ci='sd',kind='line')

# override the legend with a nicer display of the inocula
plt.gca().legend(title='Pa inoculum (cfu)',labels=['$10$','$100$','$10^3$','$10^4$']);<jupyter_output><empty_output><jupyter_text>So far, we have been using `kind='line'` to show a [`lineplot`](https://seaborn.pydata.org/generated/seaborn.lineplot.htmlseaborn.lineplot). Just like with `catplot`, you can change the `kind` of a `relplot` to to `'scatter'` to show a [`scatterplot`](https://seaborn.pydata.org/generated/seaborn.scatterplot.htmlseaborn.scatterplot) instead of a `lineplot`:<jupyter_code>from matplotlib.colors import LogNorm
sns.relplot(data=data.query("strain == 'PA103' & Ac_inoculum == 0 & time < 20"),
            x='time',y='OD600',hue='Pa_inoculum',     
            hue_norm=LogNorm(),

            aspect = 1.5,
            palette='viridis', 
            legend=False,
            ci='sd',kind='scatter')<jupyter_output><empty_output><jupyter_text>Finally, there are other aesthetics with `relplot`; you can map continuous or categorical variables to the `size` of the lines (or the size of the points for `scatterplot`); you can map categorical variables to the `style` to change the line style (or marker style for `scatterplot`):<jupyter_code>from matplotlib.colors import LogNorm
sns.relplot(data=data.query("strain == 'PA103' & Ac_inoculum in [0, 100] & time < 20"),
            x='time',y='OD600',hue='Pa_inoculum',
            
            # plot different inocula of _Ac_ with different line styles
            style='Ac_inoculum',
            hue_norm=LogNorm(),

            aspect = 1.5,
            palette='viridis', 
            legend='full',
            ci='sd',kind='line')<jupyter_output><empty_output><jupyter_text>Remember, you can also get help on `relplot` here:<jupyter_code>sns.relplot?<jupyter_output><empty_output><jupyter_text>***Let's compare PA103 across several MOIs to PA14 FliF. Just like `catplot`, `relplot` can also make "facets," showing different values of a (categorical) variable in different plots.<jupyter_code>g = sns.relplot(data=data,col='MOI',row='strain',x='time',y='OD600',ci='sd',kind='line')<jupyter_output><empty_output><jupyter_text>**Exercise**: instead of showing the two strains in separate rows, plot them as different-colored traces on the same row of plots  *****Exercise**: Let's compare the behavior of the same strain at different MOIs, by plotting each MOI as a different-colored trace; here's some code to get you started:<jupyter_code># slight hack: convert the MOI into a string so that the MOI = infinity doesn't mess up color scaling
data_cat = data.copy().sort_values(by='MOI')
data_cat.MOI = data_cat.MOI.astype('str') + ' Pa/Ac'

g = sns.relplot(data=data_cat,
                # ??
                # ??
                x='time',
                y='OD600',
                ci='sd',
                palette='viridis',
                height=5, kind='line')<jupyter_output><empty_output><jupyter_text>Interesting! Notice that MOI = infinity (e.g. no _Ac_) is in the middle, and the bacteria grow faster with some amoeba than with no amoeba. Let's use a different color scale to emphasize this.<jupyter_code>g = sns.relplot(data=data_cat,
                col='strain',
                hue='MOI',
                x='time',
                y='OD600',
                ci='sd',
                palette='inferno_r',
                height=5, kind='line')<jupyter_output><empty_output><jupyter_text>I should mention that [color scales must be chosen carefully](https://matplotlib.org/users/colormaps.html) in order to avoid distorting readers (and your own!) perception of your data, which can lead to you drawing incorrect conclusions. Let's zoom in on the region of interest: we can use the matplotlib functions [`xlim`](https://matplotlib.org/api/_as_gen/matplotlib.pyplot.xlim.html) and [`ylim`](https://matplotlib.org/api/_as_gen/matplotlib.pyplot.ylim.html) to set the limits of the `x` and `y` axes.<jupyter_code>g = sns.relplot(data=data_cat,
                col='strain',
                hue='MOI', 
                x='time',
                y='OD600',
                ci='sd',
                palette='inferno_r',
                
                # increase the aspect ratio to stretch the subplots
                aspect=1.15,
                # set the overall height to increase the size of the subplots
                height=5, 
                kind='line')

# set the limits of the x- and y-axes
plt.xlim([0,40]);
plt.ylim([0,1.]);<jupyter_output><empty_output><jupyter_text>Show different inocula of _P. aeruginosa_ as different facets, with different inocula of _A. castellanii_ as different-colored traces.<jupyter_code>g = sns.relplot(data=data_cat,
                col='strain',
                row='Pa_inoculum',
                hue='Ac_inoculum', 
                hue_norm=LogNorm(),
                
                x='time',
                y='OD600',
                ci='sd',
                palette='viridis',
                
                legend='brief',
                aspect=1.15,
                height=4, 
                kind='line')

# set a title for each subplot
g.set_titles("{col_name} / {row_name} cfu")

# set the limits of the x- and y-axes
plt.xlim([10,30]);
plt.ylim([0,0.8]);<jupyter_output><empty_output><jupyter_text>**Exercise**: plot inocula of _P. aeruginosa_ (`'Pa_inoculum'`) as colored-traces and inocula of _A. castellanii_ (`'Ac_inoculum'`) in rows of facets; plot different strains of _P. aeruginosa_ (`'strain'`) in different columns.<jupyter_code>g = sns.relplot(data=data_cat,
                # ??
                # ??
                # ??
                hue_norm=LogNorm(),
                
                x='time',
                y='OD600',
                ci='sd',
                palette='viridis',
                
                legend='brief',
                # increase the aspect ratio to stretch the subplots
                aspect=1.15,
                # set the overall height to increase the size of the subplots
                height=4, 
                kind='line')

# set the limits of the x- and y-axes
plt.xlim([10,30]);
plt.ylim([0,0.8]);<jupyter_output><empty_output><jupyter_text>A peek under the hood: `FacetGrid`Both `relplot` and `catplot` are based on a Seaborn function called [`FacetGrid`](https://seaborn.pydata.org/generated/seaborn.FacetGrid.html). You can recreate the behavior of `relplot` with `kind='line'` using code like this:<jupyter_code>g = sns.FacetGrid(data=data,col='MOI',row='strain')
g.map_dataframe(sns.lineplot,'time','OD600',ci='sd')<jupyter_output><empty_output><jupyter_text>Using `FacetGrid` is a little different than using `catplot`, because `FacetGrid` is more general than `catplot`---rather than just categorical plots, `FacetGrid` can actually be used for any type of plot! Because of this, using `FacetGrid` requires two steps:1. Create the FacetGrid, and specify how to apply variables to `row`s and `col`umns of the grid, as well as what variable(s) to use to determine the `hue`    ```python    g = sns.FacetGrid(data=data,col='MOI',row='strain')    ```2. Apply (`map`) some plotting function to each position in the grid; specify the variables to pass to this plotting function here---generally the two variables will correspond to `x` and `y`.    ```python    g.map_dataframe(sns.lineplot,'time','OD600',ci='sd')    ```        Notice that you do **not** write `x='time'`, `y='OD600'` like you did for `catplot`.  One place this can be very useful is for overlaying different types of plots; you can just call `map_dataframe` multiple times for the different styles of plots. For example, let's display a spaghetti plot of different traces on top of a mean and confidence interval:<jupyter_code>g = sns.FacetGrid(data=data.query('time < 75'),col='MOI',row='strain')
# plot a mean with a confidence interval band
g.map_dataframe(sns.lineplot,'time','OD600',ci='sd',linewidth=2,color='purple')
# show a spaghetti plot with traces of each individual well
g.map_dataframe(sns.lineplot,'time','OD600',estimator=None,units='well',alpha=0.7,linewidth=0.5)<jupyter_output><empty_output><jupyter_text>Case 3: mutants + TungstateGoal: this experiment compares mutants of PAO1 grown in the presence or absence of tungstateData: single time series.Data wrangling is very similar to what we've seen before<jupyter_code># load data from Excel file
data = pd.read_excel('./data/2018-03-29_CG007e.xlsx',
              sheet_name=0, header=36, skipfooter=4, index_col=0)

# Extract temperature vs. time as separate variable, then drop (delete)
# it from the main data frame
temps = data.loc['Temp. [°C]',:]
data = data.drop('Temp. [°C]')

# collect data into "tidy" format, with one row per observation, 
# one column for each variable (well, time, OD600)
data.columns = data.columns.rename('time')
data.index = data.index.rename('well')
data = data.reset_index()
data = data.melt(id_vars=['well'],value_name='OD600')

# build platemap with relevant metadata
platemap = plates.spec96to384(plates.prog2spec({
    "A:A,H:H,1:1,12:12": { "sterile":1 },
    "B2:G11":            { "sterile":0 },
    
    "B2:G3":   { "strain":"PAO1" , "phenotype":"susceptible" },
    "B4:G5":   { "strain":"30"   , "phenotype":"susceptible" },
    "B6:G7":   { "strain":"31"   , "phenotype":"susceptible" },
    "B8:G9":   { "strain":"35"   , "phenotype":"susceptible" },
    
    "2:2,4:4,6:6,8:8,10:10": { "tungstate":0  },
    "3:3,5:5,7:7,9:9,11:11": { "tungstate":10 },
    
    "B:B": { "dilution":2**-0 },
    "C:C": { "dilution":2**-1 },
    "D:D": { "dilution":2**-2 },
    "E:E": { "dilution":2**-3 },
    "F:F": { "dilution":2**-4 },
    "G:G": { "dilution":2**-5 }
}, include_row_column=True))

# join platemap to the raw data
data = pd.merge(data,platemap, how='inner', left_on='well',right_index=True)


# convert time series data from seconds to hours
data['time'] = data['time']/3600.0

# anything not specifically marked as sterile was innoculated
# data = data.fillna(value={'sterile':0})
data.dropna().head()
data.describe()
data.describe(exclude=[np.number])<jupyter_output><empty_output><jupyter_text>Case 3: PlottingI mostly included this to highlight an alternative plotting library, called [plotnine](http://plotnine.readthedocs.io/en/stable/index.html). Plotnine is a port of [`ggplot2`](http://ggplot2.tidyverse.org), a very popular package for R based on the ["Grammar of Graphics"](https://www.springer.com/us/book/9780387245447). This concept is very appealing, because it's much more flexible than even what we've been able to do with Seaborn.The basic idea is that a graphic is composed of **layers** stacked to create a complex plot. Each layer can map several "aesthetics" (like `x`, `y`, and `hue` as we've seen, but also `size`, `stroke`, `marker`, etc. depending on the type of layer) to different variables in the data. This allows, in principle, an amazing variety of plots to be created with remarkable simplicity. See examples [here](http://r4ds.had.co.nz/data-visualisation.html) and [here](http://r4ds.had.co.nz/graphics-for-communication.html). Pythonistas have long waited for a worthy replacement for `ggplot2`, and `plotnine` has only appeared relatively recently. I'll admit, I haven't quite gotten the hang of it—the documentation is somewhat sparse compared to Seaborn, assuming that most people using it will know how to use `ggplot2`. It seems a bit rough around the edges in a few places (it causes pandas to shout tons of warnings), but it will be interesting to watch.<jupyter_code># Seaborn plot to compare
g = sns.FacetGrid(data=data.query('sterile == 0'),
                  col='strain',row='tungstate',hue='dilution',
                  palette='viridis',legend_out=True)
g.map_dataframe(sns.lineplot,'time','OD600',legend='full',ci='sd')
g.add_legend()
from plotnine import *
import warnings
warnings.filterwarnings(action='ignore')
# numeric data needs to be cast as such, or plotnine will complain and treat it like categorical data
data = data.apply(pd.to_numeric,errors='ignore')
(ggplot(data.query('sterile == 0').dropna(), aes('time', 'OD600',color='factor(dilution)'))
 + geom_point(alpha=0.1,size=0.1)
 + scale_x_continuous(breaks=range(0,25,5))
 + scale_color_brewer()
 + stat_summary(fun_data='mean_sdl')
 + facet_grid('tungstate ~ strain'))
(ggplot(data.query('sterile == 0').dropna(), aes('time', 'OD600',color='factor(dilution)'))
 + geom_point()
 + scale_x_continuous(breaks=range(0,25,5))
 + scale_color_brewer()
 + facet_grid('tungstate ~ strain'))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12726644442637614, 'fraction_numerical': 0.025998048621839173, 'mean_word_length': 3.523862068965517, 'pattern_counts': {'":': 95, '<': 178, '<?xml version=': 0, '>': 175, 'https://': 28, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 54, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '43', 'dir': 'jupyter-structured-clean-dedup', 'id': '613068', 'license': 'CNRI-Python', 'n_tokens_mistral': 16420, 'n_tokens_neox': 15036, 'n_words': 5366, 'path': '2018-05-25_Kaz_Lab_Workshop.ipynb', 'repo_name': 'caseygrun/kazlab-python-workshop', 'size': '49196'}