{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Linear regression of (estimated) fish age to selected log-ratios for gill samples\n",
"## What does this notebook do?\n",
"This produces a scatterplot with a line showing this regression. This regression is computed twice, for two different log-ratios. These scatterplots are in the Qurro paper as Figs. 1(d) and 2(d).\n",
"\n",
"These plots, and the visualized linear regression, are only computed for gill samples: in particular, only gill samples that 1) are already used in the Qurro visualization and 2) have a valid log-ratio (i.e. do not have 0 in the numerator and/or denominator of the selected log-ratio).\n",
"\n",
"## What's with the functions?\n",
"I've moved most of the logic here to three separate functions (`load_log_ratios()`, `make_scatterplot_of_log_ratios()`, and `perform_ols_regression_on_log_ratios()`) so that replicating this sort of thing with multiple log-ratio files should be relatively painless. All you need is the TSV file with sample plot data exported from Qurro (these are in this directory as `logratios.tsv` and `10_logratios.tsv` for the `Shewanella : Synechococcales` and `Shewanella : Bottom ~10%` log-ratios, respectively).\n",
"\n",
"The reason this notebook's functionality is spread across three different functions (and not just a single one) is so that the output of each step can be easily displayed (e.g. the `head()` of a DataFrame or a drawing of a plot), since putting `.head()` in the middle of a Jupyter Notebook cell doesn't actually display the head of the DataFrame. There's almost certainly a better way to structure this code, but I figure this should be a reasonable enough solution.\n",
"\n",
"## Why do you write this documentation in a Q&A format? Seriously, reading this is weird.\n",
"Wow! That sounds like a loaded question, so I'm not going to answer it."
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {},
"outputs": [],
"source": [
"import pandas as pd\n",
"from sklearn.linear_model import LinearRegression\n",
"from matplotlib import pyplot\n",
"\n",
"def load_log_ratios(tsv_filename):\n",
" \"\"\"Returns a DataFrame of sample plot data subset to gill samples with a valid log-ratio.\"\"\"\n",
" sp = pd.read_csv(tsv_filename, sep='\\t', index_col=0, na_filter=False)\n",
" print(\"{} samples total in the visualization.\".format(len(sp.index)))\n",
" \n",
" # Remove samples with an invalid log-ratio\n",
" sp_valid = sp.loc[sp[\"Current_Natural_Log_Ratio\"] != \"null\"]\n",
" print(\"{} samples in the visualization have a valid selected log-ratio.\".format(\n",
" len(sp_valid.index)\n",
" ))\n",
" \n",
" # Remove samples that aren't from fish gills\n",
" sp_valid_gill = sp_valid.loc[sp_valid[\"sample_type_body_site\"] == \"fish gill\"]\n",
" print(\"{} samples in the visualization with a valid selected log-ratio are gill samples.\".format(\n",
" len(sp_valid_gill.index)\n",
" ))\n",
" \n",
" # We need to manually set the age and log-ratio columns as numeric for plotting them. Commenting these\n",
" # lines out will cause things to break when plotting these values.\n",
" # (Setting the columns of the DataFrame in this way results in pandas giving you a\n",
" # SettingWithCopyWarning, but after a few minutes of messing around with this I haven't figured out a\n",
" # way to get around that...)\n",
" sp_valid_gill[\"age_2\"] = pd.to_numeric(sp_valid_gill[\"age_2\"])\n",
" sp_valid_gill[\"Current_Natural_Log_Ratio\"] = pd.to_numeric(sp_valid_gill[\"Current_Natural_Log_Ratio\"])\n",
" \n",
" return sp_valid_gill\n",
"\n",
"def make_scatterplot_of_log_ratios(sp_valid_gill):\n",
" \"\"\"Draws a basic scatter-plot of estimated age vs. log-ratio for just gill samples.\n",
" \n",
" Useful as a basic sanity check before we compute the linear regression.\n",
" \"\"\"\n",
" \n",
" # We use color and transparency values that match what's used in this scatterplot in Qurro\n",
" # (this color is one of the defaults in the tableau10 colorscheme --\n",
" # see https://vega.github.io/vega/docs/schemes/#tableau10)\n",
" sp_valid_gill.plot(kind=\"scatter\", x=\"age_2\", y=\"Current_Natural_Log_Ratio\", c=\"#e45756\", alpha=0.7)\n",
" \n",
"def perform_ols_regression_on_log_ratios(sp_valid_gill, output_figure_filename):\n",
" \"\"\"Draws a scatter-plot of estimated age vs. log-ratio, with a line showing a linear regression.\n",
" \n",
" The linear regression is computed using scikit-learn. Our use of LinearRegression here is based on\n",
" this example: https://towardsdatascience.com/linear-regression-in-6-lines-of-python-5e1d0cd05b8d.\n",
" \"\"\"\n",
" # Reshape the data to get into a format that LinearRegression accepts\n",
" x = sp_valid_gill[\"age_2\"].values.reshape(-1, 1)\n",
" y = sp_valid_gill[\"Current_Natural_Log_Ratio\"].values.reshape(-1, 1)\n",
"\n",
" num_samples = len(sp_valid_gill.index)\n",
" \n",
" # Perform linear regression\n",
" lr = LinearRegression().fit(x, y)\n",
" lr_y = lr.predict(x)\n",
" \n",
" # Compute the R^2 value of this regression model\n",
" r2 = lr.score(x, y)\n",
" print(\"R^2 = {}\".format(r2))\n",
"\n",
" # Plot the same scatterplot as before, but visualize the linear regression alongside the samples\n",
" # NOTE: i've been having trouble getting transparency working, so these might end up being opaque in\n",
" # the final published figures :(\n",
" sp_valid_gill.plot(kind=\"scatter\", x=\"age_2\", y=\"Current_Natural_Log_Ratio\", c=\"#e45756\", alpha=0.7)\n",
" pyplot.plot(x, lr_y, color=\"black\")\n",
"\n",
" # Make the titles informative (and large)\n",
" pyplot.xlabel(\"Estimated fish age (age_2)\", fontsize=\"x-large\")\n",
" pyplot.ylabel(\"Current Natural Log-Ratio\", fontsize=\"x-large\")\n",
" # The double-backslashes in \\\\approx are needed in order to prevent the \\a from being interpreted\n",
" # as an escape character. Using a raw string (r\"...\") would also circumvent this problem; shout out\n",
" # to https://stackoverflow.com/a/24739481/10730311 for help in getting around this.\n",
" pyplot.title(\"Gill Samples with a Valid Log-Ratio ($n = {}$)\\n$R^2 \\\\approx {:.4f}$\".format(num_samples, r2), fontsize=\"x-large\")\n",
"\n",
" # Make ticks consistent with the other plots from Qurro (edited programmatically in Vega-Lite using the\n",
" # biggify.py script included in this repository)\n",
" pyplot.yticks([-10, -5, 0, 5])\n",
" pyplot.xticks([0, 1, 2, 3, 4])\n",
"\n",
" # Also, show gridlines (like in Fig. 1(c))\n",
" pyplot.grid(True)\n",
"\n",
" # Export the plot\n",
" pyplot.savefig(output_figure_filename, dpi=600)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Log-ratio of `Shewanella` to `Synechococcales`"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"639 samples total in the visualization.\n",
"285 samples in the visualization have a valid selected log-ratio.\n",
"143 samples in the visualization with a valid selected log-ratio are gill samples.\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"/home/marcus/Software/miniconda2/envs/skbio/lib/python3.6/site-packages/ipykernel_launcher.py:27: SettingWithCopyWarning: \n",
"A value is trying to be set on a copy of a slice from a DataFrame.\n",
"Try using .loc[row_indexer,col_indexer] = value instead\n",
"\n",
"See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy\n",
"/home/marcus/Software/miniconda2/envs/skbio/lib/python3.6/site-packages/ipykernel_launcher.py:28: SettingWithCopyWarning: \n",
"A value is trying to be set on a copy of a slice from a DataFrame.\n",
"Try using .loc[row_indexer,col_indexer] = value instead\n",
"\n",
"See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy\n"
]
},
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" Current_Natural_Log_Ratio | \n",
" age_2 | \n",
" sample_type_body_site | \n",
"
\n",
" \n",
" | Sample ID | \n",
" | \n",
" | \n",
" | \n",
"
\n",
" \n",
" \n",
" \n",
" | 11721.s38.19.gill | \n",
" 1.386294 | \n",
" 1.822654 | \n",
" fish gill | \n",
"
\n",
" \n",
" | 11721.ep25ul.1000nl.s12.5.gill | \n",
" 1.342698 | \n",
" 1.772148 | \n",
" fish gill | \n",
"
\n",
" \n",
" | 11721.s38.13.gill | \n",
" 6.555072 | \n",
" 0.864577 | \n",
" fish gill | \n",
"
\n",
" \n",
" | 11721.s38.16.gill | \n",
" 6.540870 | \n",
" 2.304387 | \n",
" fish gill | \n",
"
\n",
" \n",
" | 11721.echo5ul.200nl.s9.6.gill | \n",
" 6.877222 | \n",
" 1.772148 | \n",
" fish gill | \n",
"
\n",
" \n",
"
\n",
"
\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" Current_Natural_Log_Ratio | \n",
" age_2 | \n",
" sample_type_body_site | \n",
"
\n",
" \n",
" | Sample ID | \n",
" | \n",
" | \n",
" | \n",
"
\n",
" \n",
" \n",
" \n",
" | 11721.s38.19.gill | \n",
" 2.858111 | \n",
" 1.822654 | \n",
" fish gill | \n",
"
\n",
" \n",
" | 11721.ep25ul.1000nl.s12.5.gill | \n",
" 6.637258 | \n",
" 1.772148 | \n",
" fish gill | \n",
"
\n",
" \n",
" | 11721.s38.13.gill | \n",
" 7.471363 | \n",
" 0.864577 | \n",
" fish gill | \n",
"
\n",
" \n",
" | 11721.s38.16.gill | \n",
" 8.044947 | \n",
" 2.304387 | \n",
" fish gill | \n",
"
\n",
" \n",
" | 11721.echo5ul.200nl.s9.6.gill | \n",
" 6.877222 | \n",
" 1.772148 | \n",
" fish gill | \n",
"
\n",
" \n",
"
\n",
"