03-download-public-data.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Download public data and perform simple data analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. Google GSE41265, which will show you a link to GEO database\n",
    "\n",
    "2. Or follow this [link](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE41265) to jump directly to the GEO page for this data. Scroll down to the bottom in supplemental material. And download the link for the table called GSE41265_allGenesTPM.txt.gz.\n",
    "\n",
    "3. We also need the link to the metadata. It is <a href=\"ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE41nnn/GSE41265/matrix/\"> Series Matrix</a>. Download the file called GSE41265_series_matrix.txt.gz. \n",
    "\n",
    "4. The data is described in this reference\n",
    "* [Single-cell transcriptomics reveals bimodality in expression and splicing in immune cells](http://www.ncbi.nlm.nih.gov/pubmed/23685454) Shalek and Satija, *et al*.,  *Nature*, 2013"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reading in the data file\n",
    "\n",
    "To read the gene expression matrix, we'll use \"`pandas`\" a Python package for \"Panel Data Analysis\" (as in panels of data), which is a fantastic library for working with dataframes, and is Python's answer to R's dataframes. We'll take this opportunity to import ALL of the python libaries that we'll use today.\n",
    "\n",
    "We'll be using several additional libraries in Python:\n",
    "\n",
    "3. [`matplotlib`](http://matplotlib.org/) - This is the base plotting library in Python.\n",
    "1. [`numpy`](http://www.numpy.org/) - (pronounced \"num-pie\") which is basis for most scientific packages. It's basically a nice-looking Python interface to C code. It's very fast.\n",
    "2. [`pandas`](http://pandas.pydata.org) - This is the \"DataFrames in Python.\" (like R's nice dataframes) They're a super convenient form that's based on `numpy` so they're fast. And you can do convenient things like calculate mea n and variance very easily.\n",
    "4. [`scipy`](http://www.scipy.org/) - (pronounced \"sigh-pie\") \"Scientific Python\" - Contains statistical methods and calculations\n",
    "4. [`seaborn`](http://web.stanford.edu/~mwaskom/software/seaborn/index.html) - Statistical plotting library. To be completely honest, R's plotting and graphics capabilities are much better than Python's. However, Python is a really nice langauge to learn and use, it's very memory efficient, can be parallized well, and has a very robust machine learning library, `scikit-learn`, which has a very nice and consistent interface. So this is Python's answer to `ggplot2` (very popular R library for plotting) to try and make plotting in Python nicer looking and to make statistical plots easier to do.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Python plotting library\n",
    "import matplotlib.pyplot as plt\n",
    "# This is necessary to show the plotted figures inside the notebook -- \"inline\" with the notebook cells\n",
    "%matplotlib inline\n",
    "\n",
    "# Numerical python library (pronounced \"num-pie\")\n",
    "import numpy as np\n",
    "\n",
    "# Dataframes in Python\n",
    "import pandas as pd\n",
    "\n",
    "# Statistical plotting library we'll use\n",
    "import seaborn as sns"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll read in the data using `pandas` and look at the first 5 rows of the dataframe with the dataframe-specific function `.head()`. Whenever I read a new table or modify a dataframe, I **ALWAYS** look at it to make sure it was correctly imported and read in, and I want you to get into the same habit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read the data table\n",
    "# You may need to change the path to the file (what's in quotes below) relative \n",
    "# to where you downloaded the file and where this notebook is\n",
    "shalek2013_expression = pd.read_table('data/GSE41265_allGenesTPM.txt.gz', \n",
    "                               \n",
    "                                     # Sets the first (Python starts counting from 0 not 1) column as the row names\n",
    "                                      index_col=0, \n",
    "\n",
    "                                     # Tells pandas to decompress the gzipped file\n",
    "                                      compression='gzip')\n",
    "\n",
    "\n",
    "print(shalek2013_expression.shape)\n",
    "shalek2013_expression.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So we have 21 columns but looks like here `pandas` by default is showing a maximum of 20 so let's change the setting so we can see ALL of the samples instead of just skipping single cell 11 (**S11**). Let's change to 50 for good measure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pd.options.display.max_columns = 50\n",
    "pd.options.display.max_rows = 50\n",
    "shalek2013_expression.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can see all the samples!\n",
    "\n",
    "Let's take a look at the full size of the matrix with `.shape`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "shalek2013_expression.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reading in the metadata"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "shalek2013_metadata = pd.read_table('data/GSE41265_series_matrix.txt.gz',\n",
    "                                    compression = 'gzip',\n",
    "                                    skiprows=33, \n",
    "                                    index_col=0)\n",
    "print(shalek2013_metadata.shape)\n",
    "shalek2013_metadata"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's transpose this matrix so the samples are the rows, and the features are the columns. We'll do that with `.T`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "shalek2013_metadata = shalek2013_metadata.T\n",
    "shalek2013_metadata"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we'll do some mild data cleaning. Notice that the columns have the exclamation point at the beginning, so let's get rid of that. In computer science, you keep letters between quotes, and you call those \"strings.\" Let's talk about the string function `.strip()`. This removes any characters that are on the outer edges of the string. For example, let's take the string \"Whoooo!!!!!!!\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can access the column names with `dataframe.columns`, like below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "shalek2013_metadata.columns"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can map the stripping function to every item of the columns. In Python, the square brackets (`[` and `]`) show that we're making a list. What we're doing below is called a \"list comprehension.\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "[x.strip('!') for x in shalek2013_metadata.columns]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In `pandas`, we can do the same thing by `map`-ping a `lambda`, which is a small, anonymous function that does one thing. It's called \"anonymous\" because it doesn't have a name. `map` runs the function on every element of the columns."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "shalek2013_metadata.columns.map(lambda x: x.strip('!'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above `lambda` is the same as if we had written a *named* function called `remove_exclamation`, as below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def remove_exclamation(x):\n",
    "    return x.strip('!')\n",
    "\n",
    "shalek2013_metadata.columns.map(remove_exclamation)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can assign the new column names to our matrix:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "shalek2013_metadata.columns = shalek2013_metadata.columns.map(lambda x: x.strip('!'))\n",
    "shalek2013_metadata.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Okay, now we're ready to do some analysis! \n",
    "\n",
    "We've looked at the top of the dataframe by using `head()`. By default, this shows the first 5 rows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbgrader": {
     "checksum": "890994ad468ad7529d7baa0a4bbd2dbf",
     "grade": false,
     "grade_id": "5_head_df",
     "locked": false,
     "solution": false
    }
   },
   "outputs": [],
   "source": [
    "shalek2013_expression.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "checksum": "ef2e6d24e06d536e932949e2cbd40f75",
     "grade": false,
     "grade_id": "5_head_df_n",
     "locked": false,
     "solution": false
    }
   },
   "source": [
    "To specify a certain number of rows, put a number between the parentheses."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbgrader": {
     "checksum": "a600224e6ea5acb0981cbb7cc6ab62d1",
     "grade": false,
     "grade_id": "5_head_df_8",
     "locked": false,
     "solution": false
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "shalek2013_expression.head(8)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "checksum": "b8ac47c7ac6c55e00f86257e4f42ffc8",
     "grade": false,
     "grade_id": "explain_seaborn_boxplot",
     "locked": false,
     "solution": false
    }
   },
   "source": [
    "Let's get a sense of this data by plotting the distributions using `boxplot` from seaborn. To save the output, we'll need to get access to the current figure, and save this to a variable using `plt.gcf()`. And then we'll save this figure with `fig.savefig(\"filename.pdf\")`. You can use other extensions (e.g. \"`.png`\", \"`.tiff`\" and it'll automatically save as that forma)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbgrader": {
     "checksum": "32fbdfc1bce49aa639d3c85d2e6bc2c9",
     "grade": false,
     "grade_id": "run_seaborn_boxplot",
     "locked": false,
     "solution": false
    }
   },
   "outputs": [],
   "source": [
    "sns.boxplot(shalek2013_expression['S1'])\n",
    "\n",
    "# gcf = Get current figure\n",
    "fig = plt.gcf()\n",
    "fig.savefig('shalek2013_expression_boxplot.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "checksum": "1e356f3c2f978b3564d78054d371f457",
     "grade": false,
     "grade_id": "explain_log_scale",
     "locked": false,
     "solution": false
    }
   },
   "source": [
    "Notice the 140,000 maximum ... Oh right we have expression data and the scales are enormous... Let's add 1 to all values and take the log2 of the data. We add one because log(0) is undefined and then all our logged values start from zero too. This \"$\\log_2(TPM + 1)$\" is a very common transformation of expression data so it's easier to analyze."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbgrader": {
     "checksum": "3c385e7b0e0d13a047438fb16ab575cf",
     "grade": false,
     "grade_id": "do_log_scale",
     "locked": false,
     "solution": false
    }
   },
   "outputs": [],
   "source": [
    "expression_logged = np.log2(shalek2013_expression+1)\n",
    "expression_logged.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbgrader": {
     "checksum": "0591be7ce718355d190d114034f18a1e",
     "grade": false,
     "grade_id": "seaborn_logged_expression",
     "locked": false,
     "solution": false
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "sns.boxplot(x=expression_logged['S1'])\n",
    "\n",
    "# gcf = Get current figure\n",
    "#fig = plt.gcf()\n",
    "#fig.savefig('expression_logged_boxplot.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "checksum": "db966438110ed2dae2cdf6258e031549",
     "grade": false,
     "grade_id": "filtering",
     "locked": false,
     "solution": false
    }
   },
   "source": [
    "## Filtering expression data\n",
    "\n",
    "Seems like a lot of genes are near zero, which means we need to filter our genes.\n",
    "\n",
    "We can ask which genes have log2 expression values are less than 2 (weird example I know - stay with me). This creates a dataframe of `boolean` values of True/False."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbgrader": {
     "checksum": "d7bb3501678238f1a41465f9e16659ac",
     "grade": false,
     "grade_id": "boolean_matrix",
     "locked": false,
     "solution": false
    }
   },
   "outputs": [],
   "source": [
    "at_most_2 = expression_logged < 2\n",
    "at_most_2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "checksum": "fbfc82a628d9e94a5cd9d1e782bc0411",
     "grade": false,
     "grade_id": "explain_boolean_df",
     "locked": false,
     "solution": false
    }
   },
   "source": [
    "What's nice about booleans is that False is 0 and True is 1, so we can sum to get the number of \"Trues.\" This is a simple, clever way that we can filter on a count for the data. We **could** use this boolean dataframe to filter our original dataframe, but then we lose information. For all values that are greater than 2, it puts in a \"not a number\" - \"NaN.\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbgrader": {
     "checksum": "32481371456e30b37ce803a2acd73ab1",
     "grade": false,
     "grade_id": "filter_with_boolean_df",
     "locked": false,
     "solution": false
    }
   },
   "outputs": [],
   "source": [
    "expression_at_most_2 = expression_logged[expression_logged < 2]\n",
    "print(expression_at_most_2.shape)\n",
    "expression_at_most_2.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "checksum": "cda74abb93b8b4f13aba3d463b3a406b",
     "grade": false,
     "grade_id": "ex3_question",
     "locked": false,
     "solution": false
    }
   },
   "source": [
    "### Exercise 4: Crude filtering on expression data\n",
    "\n",
    "Create a dataframe called \"`expression_greater_than_5`\" which contains only values that are greater than 5 from `expression_logged`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "expression_logged.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "expression_greater_than_5 = expression_logged[expression_logged > 5]\n",
    "expression_greater_than_5.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "checksum": "2db5edb4cd9ecc42d78f1a8053a233fb",
     "grade": false,
     "grade_id": "smarter_filtering",
     "locked": false,
     "solution": false
    }
   },
   "source": [
    "\n",
    "The crude filtering above is okay, but we're smarter than that. We want to use the filtering in the paper: \n",
    "\n",
    "> *... discarded genes that were not appreciably expressed (transcripts per million (TPM) > 1) in at least three individual cells, retaining 6,313 genes for further analysis.*\n",
    "\n",
    "We want to do THAT, but first we need a couple more concepts. The first one is summing booleans.\n",
    "\n",
    "## A smarter way to filter\n",
    "\n",
    "Remember that booleans are really 0s (`False`) and 1s (`True`)? This turns out to be VERY convenient and we can use this concept in clever ways.\n",
    "\n",
    "We can use `.sum()` on a boolean matrix to get the number of genes with expression greater than 10 for each sample:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbgrader": {
     "checksum": "f767b67c376b04e51ce27aa5e321c7ef",
     "grade": false,
     "grade_id": "boolean_sum",
     "locked": false,
     "solution": false
    }
   },
   "outputs": [],
   "source": [
    "(expression_logged > 10).sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "checksum": "f92f2a650571923f467ec0421c3ff87c",
     "grade": false,
     "grade_id": "explain_sum_axis1",
     "locked": false,
     "solution": false
    }
   },
   "source": [
    "`pandas` is column-oriented and by default, it will give you a sum for each column. But **we** want a sum for each row. How do we do that?\n",
    "\n",
    "\n",
    "We can sum the boolean matrix we created with \"`expression_logged < 10`\" along `axis=1` (along the samples) to get **for each gene, how many samples have expression less than 10**. In `pandas`, this column is called a \"`Series`\" because it has only one dimension - its length. Internally, `pandas` stores dataframes as a bunch of columns - specifically these `Series`ssssss.\n",
    "\n",
    "This turns out to be not that many."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbgrader": {
     "checksum": "9b32cd77e232e4d0c37f285bc2b57d32",
     "grade": false,
     "grade_id": "boolean_sum_axis1",
     "locked": false,
     "solution": false
    }
   },
   "outputs": [],
   "source": [
    "(expression_logged > 10).sum(axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "checksum": "ac97d0d811684014dba3a24c502b5c4e",
     "grade": false,
     "grade_id": "ex4_question",
     "locked": false,
     "solution": false
    }
   },
   "source": [
    "Now we can apply ANOTHER filter and find genes that are \"present\" (expression greater than 10) in at least 5 samples. We'll save this as the variable `genes_of_interest`. Notice that this doesn't the `genes_of_interest` but rather the list at the bottom. This is because what you see under a code cell is the output of the last thing you called. The \"hash mark\"/\"number sign\" \"`#`\" is called a **comment character** and makes the rest of the line after it not read by the Python language.\n",
    "\n",
    "### Exercise 5: Commenting and uncommenting\n",
    "\n",
    "To see `genes_of_interest`, \"uncomment\" the line by removing the hash sign, and commenting out the list `[1, 2, 3]`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbgrader": {
     "checksum": "00ecfb44b16b5d9c3b57f6d4c9bb8913",
     "grade": false,
     "grade_id": "ex4_answer",
     "locked": false,
     "solution": false
    }
   },
   "outputs": [],
   "source": [
    "genes_of_interest = (expression_logged > 10).sum(axis=1) >= 5\n",
    "#genes_of_interest\n",
    "[1, 2, 3]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "checksum": "9c4e6f3fc04066964e983b5c0652d943",
     "grade": false,
     "grade_id": "get_rows",
     "locked": false,
     "solution": false
    }
   },
   "source": [
    "## Getting only rows that you want (aka subsetting)\n",
    "\n",
    "Now we have some genes that we want to use - how do you pick just those? This can also be called \"subsetting\" and in `pandas` has the technical name [indexing](http://pandas.pydata.org/pandas-docs/stable/indexing.html)\n",
    "\n",
    "In `pandas`, to get the rows (genes) you want using their name (gene symbol) or boolean matrix, you use `.loc[rows_you_want]`. Check it out below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbgrader": {
     "checksum": "6ee77eb36ce8d6d751aef7b61d3401a1",
     "grade": false,
     "grade_id": "show_row_Subsetting",
     "locked": false,
     "solution": false
    }
   },
   "outputs": [],
   "source": [
    "expression_filtered = expression_logged.loc[genes_of_interest]\n",
    "print(expression_filtered.shape)  # shows (nrows, ncols) - like in manhattan you do the Street then the Avenue\n",
    "expression_filtered.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "checksum": "111d70a099561c8748ddd1aa98659f58",
     "grade": false,
     "grade_id": "ex5_question",
     "locked": false,
     "solution": false
    }
   },
   "source": [
    "Wow, our matrix is very small - 197 genes! We probably don't want to filter THAT much... I'd say a range of 5,000-15,000 genes after filtering is a good ballpark. Not too big so it's impossible to work with but not too small that you can't do any statistics.\n",
    "\n",
    "We'll get closer to the expression data created by the paper. Remember that they filtered on genes that had expression greater than 1 in at least 3 *single cells*. We'll filter for expression greater than 1 in at least 3 *samples* for now - we'll get to the single stuff in a bit. For now, we'll filter on all samples.\n",
    "\n",
    "### Exercise: Filtering on the presence of genes\n",
    "\n",
    "Create a dataframe called `expression_filtered_by_all_samples` that consists only of genes that have expression greater than 1 in at least 3 samples.\n",
    "\n",
    "#### Hint for `IndexingError: Unalignable boolean Series key provided`\n",
    "\n",
    "If you're getting this error, double-check your `.sum()` command. Did you remember to specify that you want to get the number of cells (columns) that express each **gene** (row)? Remember that `.sum()` by default gives you the sum over columns, but since genes are the rows .... How do you get the sum over rows?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "genes_of_interest = (expression_logged > 1).sum(axis=1) >= 3\n",
    "\n",
    "expression_filtered_by_all_samples = expression_logged.loc[genes_of_interest]\n",
    "print(expression_filtered_by_all_samples.shape)\n",
    "expression_filtered_by_all_samples.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "checksum": "a6422e612e68ed72b06e5d6e63eee909",
     "grade": false,
     "grade_id": "explain_boxplot_filtered_all_samples",
     "locked": false,
     "solution": false
    }
   },
   "source": [
    "Just for fun, let's see how our the distributions in our expression matrix have changed. If you want to save the figure, you can:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbgrader": {
     "checksum": "0760115865abdd218d1695de1c5c2ca4",
     "grade": false,
     "grade_id": "do_boxplot_filtered_all_samples",
     "locked": false,
     "solution": false
    }
   },
   "outputs": [],
   "source": [
    "sns.boxplot(expression_filtered_by_all_samples)\n",
    "\n",
    "# gcf = Get current figure\n",
    "fig = plt.gcf()\n",
    "fig.savefig('expression_filtered_by_all_samples_boxplot.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "checksum": "0aab90d4a140c7ca74026c8b9dec50e8",
     "grade": false,
     "grade_id": "subsetting_columns",
     "locked": false,
     "solution": false
    }
   },
   "source": [
    "## Getting only the columns you want\n",
    "\n",
    "In the next exercise, we'll get just the single cells\n",
    "\n",
    "For the next step, we're going to pull out just the pooled - which are conveniently labeled as \"P#\". We'll do this using a [list comprehension](http://www.pythonforbeginners.com/basics/list-comprehensions-in-python), which means we'll create a new list based on the items in `shalek2013_expression.columns` and whether or not they start with the letter `'P'`.\n",
    "\n",
    "In Python, things in square brackets (`[]`) are lists unless indicated otherwise. We are using a list comprehension here instead of a map, because we only want a *subset* of the columns, rather than *all* of them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbgrader": {
     "checksum": "f59d52c19044bc8ba61b209a1166583f",
     "grade": false,
     "grade_id": "list_comprehension_example",
     "locked": false,
     "solution": false
    }
   },
   "outputs": [],
   "source": [
    "pooled_ids = [x for x in expression_logged.columns if x.startswith('P')]\n",
    "pooled_ids"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "checksum": "1ca93f9bd3b8abb59e49313d07360ca1",
     "grade": false,
     "grade_id": "subset_columns",
     "locked": false,
     "solution": false
    }
   },
   "source": [
    "We'll access the columns we want using this bracket notation (note that this only works for columns, not rows)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbgrader": {
     "checksum": "91d8cc680c7a1542c21a9e86c59b87ca",
     "grade": false,
     "grade_id": "show_subset_columns",
     "locked": false,
     "solution": false
    }
   },
   "outputs": [],
   "source": [
    "pooled = expression_logged[pooled_ids]\n",
    "pooled.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "checksum": "22eee373e463d204f19ecd1f6aea429c",
     "grade": false,
     "grade_id": "subset_columns_loc",
     "locked": false,
     "solution": false
    }
   },
   "source": [
    "We could do the same thing using `.loc` but we would need to put a colon \"`:`\" in the \"rows\" section (first place) to show that we want \"all rows.\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbgrader": {
     "checksum": "56f23770caf1453460a61774a41806f1",
     "grade": false,
     "grade_id": "show_subset_columns_loc",
     "locked": false,
     "solution": false
    }
   },
   "outputs": [],
   "source": [
    "expression_logged.loc[:, pooled_ids].head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "checksum": "8156a739fa505658538c660c05aa3cc8",
     "grade": false,
     "grade_id": "ex6_question",
     "locked": false,
     "solution": false
    }
   },
   "source": [
    "### Exercise: Make a dataframe of only single samples\n",
    "\n",
    "Use list comprehensions to make a list called `single_ids` that consists only of single cells, and use that list to subset `expression_logged` and create a dataframe called `singles`. (Hint - how are the single cells ids different from the pooled ids?)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "single_ids = [x for x in expression_logged.columns if x.startswith('S')]\n",
    "singles = expression_logged[single_ids]\n",
    "print(singles.shape)\n",
    "singles.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "checksum": "416d62b817c0764841de900a39ac91d9",
     "grade": false,
     "grade_id": "ex7_question",
     "locked": false,
     "solution": false
    }
   },
   "source": [
    "## Using two different dataframes for filtering\n",
    "\n",
    "### Exercise 8: Filter the full dataframe using the singles dataframe\n",
    "\n",
    "Now we'll actually do the filtering done by the paper. Using the `singles` dataframe you just created, get the genes that have expression greater than 1 in at least 3 single cells, and use that to filter `expression_logged`. Call this dataframe `expression_filtered_by_singles`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rows = (singles > 1).sum(axis=1) > 3\n",
    "\n",
    "expression_filtered_by_singles = expression_logged.loc[rows]\n",
    "print(expression_filtered_by_singles.shape)\n",
    "expression_filtered_by_singles.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "checksum": "360185ccc7c8494f9f8739fa251ed601",
     "grade": false,
     "grade_id": "expression_filtered_by_singles_boxplot",
     "locked": false,
     "solution": false
    }
   },
   "source": [
    "Let's make a boxplot again to see how the data has changed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbgrader": {
     "checksum": "8444dc3f576bcc9944d0f61d87bfc14f",
     "grade": false,
     "grade_id": "show_expression_filtered_by_singles_boxplot",
     "locked": false,
     "solution": false
    }
   },
   "outputs": [],
   "source": [
    "sns.boxplot(expression_filtered_by_singles)\n",
    "\n",
    "fig = plt.gcf()\n",
    "fig.savefig('expression_filtered_by_singles_boxplot.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "checksum": "64b76d925e25071893a41fb076ac23e9",
     "grade": false,
     "grade_id": "why_did_this_matter",
     "locked": false,
     "solution": false
    }
   },
   "source": [
    "This is much nicer because now we don't have so many zeros and each sample has a reasonable dynamic range.\n",
    "\n",
    "## Why did this filtering even matter?\n",
    "\n",
    "You may be wondering, we did all this work to remove some zeros..... so the FPKM what? Let's take a look at how this affects the relationships between samples using `sns.jointplot` from seaborn, which will plot a correlation scatterplot. This also calculates the [Pearson correlation](https://en.wikipedia.org/wiki/Pearson_product-moment_correlation_coefficient), a linear correlation metric.\n",
    "\n",
    "Let's first do this on the unlogged data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbgrader": {
     "checksum": "a41f24e4ab1315fea1b29f6bd5b1732c",
     "grade": false,
     "grade_id": "jointplot_unlogged",
     "locked": false,
     "solution": false
    }
   },
   "outputs": [],
   "source": [
    "sns.jointplot(shalek2013_expression['S1'], shalek2013_expression['S2'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "checksum": "2eff1b6e1987e0f4d3d4229442340923",
     "grade": false,
     "grade_id": "jointplot_logged_weird",
     "locked": false,
     "solution": false
    }
   },
   "source": [
    "Pretty funky looking huh? That's why we logged it :)\n",
    "\n",
    "Now let's try this on the logged data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbgrader": {
     "checksum": "32929a1501d8e7db5b439de317c4ff86",
     "grade": false,
     "grade_id": "jointplot_expression_logged",
     "locked": false,
     "solution": false
    }
   },
   "outputs": [],
   "source": [
    "sns.jointplot(expression_logged['S1'], expression_logged['S2'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "checksum": "45e02c41ba1635fad5f65ceee9174c56",
     "grade": false,
     "grade_id": "jointplot_expression_logged_changed",
     "locked": false,
     "solution": false
    }
   },
   "source": [
    "Hmm our pearson correlation increased from 0.62 to 0.64. Why could that be?\n",
    "\n",
    "Let's look at this same plot using the filtered data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbgrader": {
     "checksum": "78df5310ff41cf06d89255009571445b",
     "grade": false,
     "grade_id": "show_expression_filtered_by_singles_jointplot",
     "locked": false,
     "solution": false
    }
   },
   "outputs": [],
   "source": [
    "sns.jointplot(expression_filtered_by_singles['S1'], expression_filtered_by_singles['S2'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.3"
  },
  "name": "_merged_merged"
 },
 "nbformat": 4,
 "nbformat_minor": 1
}