diff --git a/figure_generation_clinical_day4_afermentans_reordered.ipynb b/figure_generation_clinical_day4_afermentans_reordered.ipynb
new file mode 100644
index 0000000..6111b9f
--- /dev/null
+++ b/figure_generation_clinical_day4_afermentans_reordered.ipynb
@@ -0,0 +1,1052 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": 20,
+ "id": "dd7fc51e-b713-4b42-9cfc-3ab3ba59e7cb",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# ----------------------------------------------------------------------------\n",
+ "# SkinCom Analysis - Perform t-test & boxplots for log-ratios\n",
+ "# Author: Sherlyn Weng\n",
+ "# Version: 1.0.0\n",
+ "# Mmaintainer: Sherlyn Weng\n",
+ "# Email: y1weng@ucsd.edu\n",
+ "# This code is adapted from https://github.com/knightlab-analyses/reference-frames\n",
+ "# ----------------------------------------------------------------------------"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 21,
+ "id": "515d875f-5c34-4203-a374-67ec76fbe6ba",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import numpy as np\n",
+ "import pandas as pd\n",
+ "import seaborn as sns\n",
+ "import matplotlib.pyplot as plt\n",
+ "import matplotlib as mpl\n",
+ "from biom import load_table\n",
+ "from scipy.stats import ttest_rel, wilcoxon, ttest_ind, pearsonr, spearmanr\n",
+ "from skbio.stats.composition import multiplicative_replacement\n",
+ "import matplotlib.gridspec as gridspec\n",
+ "\n",
+ "from collections import namedtuple\n",
+ "\n",
+ "%matplotlib inline\n",
+ "\n",
+ "import warnings\n",
+ "warnings.filterwarnings('ignore')\n",
+ "\n",
+ "\"\"\"\n",
+ "Here, we have implemented an extended version of scipy's paired t-test function.\n",
+ "see https://github.com/scipy/scipy/blob/master/scipy/stats/mstats_basic.py#L1085\n",
+ "\n",
+ "This should really be pushed to scipy -- any volunteers?\n",
+ "\"\"\"\n",
+ "import numpy as np\n",
+ "from scipy.stats import t\n",
+ "from scipy.special import betainc\n",
+ "from collections import namedtuple\n",
+ "\n",
+ "\n",
+ "Ttest_relResult = namedtuple('Ttest_relResult', \n",
+ " ('statistic', 'pvalue', 'lower_CI', 'upper_CI', 'ddof'))\n",
+ "\n",
+ "\n",
+ "def ttest_rel(x, y):\n",
+ " \n",
+ " D = x - y\n",
+ " n = len(D)\n",
+ " \n",
+ " sd = np.std(D, ddof=1) \n",
+ " se = sd / np.sqrt(n)\n",
+ " mu = np.mean(D)\n",
+ " T = mu / se\n",
+ " ddof = n-1\n",
+ " ci = se * t.ppf( 0.05 / 2, ddof) \n",
+ " lower = mu - np.abs(ci)\n",
+ " upper = mu + np.abs(ci)\n",
+ " \n",
+ " prob = betainc(0.5*ddof, 0.5, ddof/(ddof + T*T))\n",
+ " return Ttest_relResult(T, prob, lower, upper, int(ddof))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 22,
+ "id": "1b2dca95-1936-4547-98ec-c7e609a219e6",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def paired_t_test(category, md, df):\n",
+ " \"\"\" Performs paired t-test on a dataframe with complementing\n",
+ " sample metadata.\n",
+ "\n",
+ " Parameters\n",
+ " ----------\n",
+ " category : str\n",
+ " Category to test in the dataframe df\n",
+ " md : pd.DataFrame\n",
+ " Sample metadata specifying groupings\n",
+ " df : pd.Dataframe\n",
+ " Counts of microbes across all samples\n",
+ " \n",
+ " Returns\n",
+ " -------\n",
+ " T : float\n",
+ " T-statistic\n",
+ " prob : float\n",
+ " p-value\n",
+ " lower : float\n",
+ " lower 95% CI bound\n",
+ " upper : float\n",
+ " upper 95% CI bound\n",
+ " ddof : int\n",
+ " degrees of freedom\n",
+ " \"\"\"\n",
+ " cats = list(md[category].value_counts().index)\n",
+ " src = md[md[category] == cats[1]].index\n",
+ " print('hi')\n",
+ " print(src)\n",
+ " dest = md[md[category] == cats[0]].index\n",
+ " return ttest_rel(df.loc[dest].values, df.loc[src].values)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "5242c6a2-b6b2-497e-aa34-225f3a9c1cab",
+ "metadata": {},
+ "source": [
+ "## Input Directory"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 23,
+ "id": "0f656bb7-fbf6-4b0e-909b-d94835646df4",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "input_dir = \"..\""
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "f9b78015-1a9a-4d8e-97d0-e7ff85a1e82c",
+ "metadata": {},
+ "source": [
+ "## Variables"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 24,
+ "id": "af3b610d-80e5-4525-b1bf-86d36a88e9cb",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "syncom_strains = [\"Corynebacterium afermentans\",\n",
+ " \"Cutibacterium acnes\",\n",
+ " \"Micrococcus luteus\",\n",
+ " \"Staphylococcus aureus\",\n",
+ " \"Staphylococcus capitis\",\n",
+ " \"Staphylococcus epidermidis\",\n",
+ " \"Staphylococcus hominis\",\n",
+ " \"Staphylococcus warneri\",\n",
+ " \"Streptococcus mitis\"]"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "1164688b-1357-4d73-bedb-66dc5b71386a",
+ "metadata": {
+ "tags": []
+ },
+ "source": [
+ "## SLES"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 26,
+ "id": "da5f89a9-1fc4-417f-82ef-271c89c4c6e7",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "hi\n",
+ "Index(['post_left_sub12', 'post_left_sub13', 'post_left_sub1',\n",
+ " 'post_left_sub2'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hi\n",
+ "Index(['post_left_sub12', 'post_left_sub13', 'post_left_sub1',\n",
+ " 'post_left_sub2'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hi\n",
+ "Index(['post_left_sub12', 'post_left_sub13', 'post_left_sub1',\n",
+ " 'post_left_sub2'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hi\n",
+ "Index(['post_left_sub12', 'post_left_sub13', 'post_left_sub1',\n",
+ " 'post_left_sub2'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hi\n",
+ "Index(['post_left_sub12', 'post_left_sub13', 'post_left_sub1',\n",
+ " 'post_left_sub2'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hi\n",
+ "Index(['post_left_sub12', 'post_left_sub13', 'post_left_sub1',\n",
+ " 'post_left_sub2'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hi\n",
+ "Index(['post_left_sub12', 'post_left_sub13', 'post_left_sub1',\n",
+ " 'post_left_sub2'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hi\n",
+ "Index(['post_left_sub12', 'post_left_sub13', 'post_left_sub1',\n",
+ " 'post_left_sub2'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hi\n",
+ "Index(['post_left_sub12', 'post_left_sub13', 'post_left_sub1',\n",
+ " 'post_left_sub2'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hi\n",
+ "Index(['post_left_sub12', 'post_left_sub13', 'post_left_sub1',\n",
+ " 'post_left_sub2'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hi\n",
+ "Index(['post_left_sub12', 'post_left_sub13', 'post_left_sub1',\n",
+ " 'post_left_sub2'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hi\n",
+ "Index(['post_left_sub12', 'post_left_sub13', 'post_left_sub1',\n",
+ " 'post_left_sub2'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hi\n",
+ "Index(['post_left_sub12', 'post_left_sub13', 'post_left_sub1',\n",
+ " 'post_left_sub2'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hi\n",
+ "Index(['post_left_sub12', 'post_left_sub13', 'post_left_sub1',\n",
+ " 'post_left_sub2'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hi\n",
+ "Index(['post_left_sub12', 'post_left_sub13', 'post_left_sub1',\n",
+ " 'post_left_sub2'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hi\n",
+ "Index(['post_left_sub12', 'post_left_sub13', 'post_left_sub1',\n",
+ " 'post_left_sub2'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hi\n",
+ "Index(['post_left_sub12', 'post_left_sub13', 'post_left_sub1',\n",
+ " 'post_left_sub2'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hi\n",
+ "Index(['post_left_sub12', 'post_left_sub13', 'post_left_sub1',\n",
+ " 'post_left_sub2'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hi\n",
+ "Index(['post_left_sub12', 'post_left_sub13', 'post_left_sub1',\n",
+ " 'post_left_sub2'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hi\n",
+ "Index(['post_left_sub12', 'post_left_sub13', 'post_left_sub1',\n",
+ " 'post_left_sub2'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hi\n",
+ "Index(['post_left_sub12', 'post_left_sub13', 'post_left_sub1',\n",
+ " 'post_left_sub2'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hi\n",
+ "Index(['post_left_sub12', 'post_left_sub13', 'post_left_sub1',\n",
+ " 'post_left_sub2'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hi\n",
+ "Index(['post_left_sub12', 'post_left_sub13', 'post_left_sub1',\n",
+ " 'post_left_sub2'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hi\n",
+ "Index(['post_left_sub12', 'post_left_sub13', 'post_left_sub1',\n",
+ " 'post_left_sub2'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hi\n",
+ "Index(['post_left_sub12', 'post_left_sub13', 'post_left_sub1',\n",
+ " 'post_left_sub2'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hi\n",
+ "Index(['post_left_sub12', 'post_left_sub13', 'post_left_sub1',\n",
+ " 'post_left_sub2'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hi\n",
+ "Index(['post_left_sub12', 'post_left_sub13', 'post_left_sub1',\n",
+ " 'post_left_sub2'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hi\n",
+ "Index(['post_left_sub12', 'post_left_sub13', 'post_left_sub1',\n",
+ " 'post_left_sub2'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hi\n",
+ "Index(['post_left_sub12', 'post_left_sub13', 'post_left_sub1',\n",
+ " 'post_left_sub2'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hi\n",
+ "Index(['post_left_sub12', 'post_left_sub13', 'post_left_sub1',\n",
+ " 'post_left_sub2'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hi\n",
+ "Index(['post_left_sub12', 'post_left_sub13', 'post_left_sub1',\n",
+ " 'post_left_sub2'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hi\n",
+ "Index(['post_left_sub12', 'post_left_sub13', 'post_left_sub1',\n",
+ " 'post_left_sub2'],\n",
+ " dtype='object', name='sampleid')\n"
+ ]
+ },
+ {
+ "data": {
+ "text/html": [
+ "\n",
+ "\n",
+ "
\n",
+ " \n",
+ " \n",
+ " | \n",
+ " 24hControl_vs_24hChemical_tstat | \n",
+ " 24hControl__vs_24hChemical_pvalue | \n",
+ " 24hControl_vs_24hChemical_upperCI | \n",
+ " 24hControl_vs_24hChemical_lowerCI | \n",
+ "
\n",
+ " \n",
+ " \n",
+ " \n",
+ " | log(C.acnes/C.afermentans) | \n",
+ " -1.006584 | \n",
+ " 0.388289 | \n",
+ " -2.011125 | \n",
+ " 1.044617 | \n",
+ "
\n",
+ " \n",
+ " | log(M.luteus/C.afermentans) | \n",
+ " -0.247151 | \n",
+ " 0.820740 | \n",
+ " -2.271485 | \n",
+ " 1.944100 | \n",
+ "
\n",
+ " \n",
+ " | log(S.aureus/C.afermentans) | \n",
+ " -0.953509 | \n",
+ " 0.410676 | \n",
+ " -0.886183 | \n",
+ " 0.477579 | \n",
+ "
\n",
+ " \n",
+ " | log(S.capitis/C.afermentans) | \n",
+ " 0.229355 | \n",
+ " 0.833340 | \n",
+ " -2.769521 | \n",
+ " 3.199716 | \n",
+ "
\n",
+ " \n",
+ " | log(S.epidermidis/C.afermentans) | \n",
+ " -0.075735 | \n",
+ " 0.944398 | \n",
+ " -2.304923 | \n",
+ " 2.197770 | \n",
+ "
\n",
+ " \n",
+ " | log(S.hominis/C.afermentans) | \n",
+ " 0.859905 | \n",
+ " 0.453072 | \n",
+ " -0.345936 | \n",
+ " 0.602096 | \n",
+ "
\n",
+ " \n",
+ " | log(S.warneri/C.afermentans) | \n",
+ " -1.075817 | \n",
+ " 0.360825 | \n",
+ " -1.542561 | \n",
+ " 0.763129 | \n",
+ "
\n",
+ " \n",
+ " | log(S.mitis/C.afermentans) | \n",
+ " 0.287471 | \n",
+ " 0.792466 | \n",
+ " -0.252952 | \n",
+ " 0.303188 | \n",
+ "
\n",
+ " \n",
+ "
\n",
+ "
"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "# load md\n",
+ "metadata = pd.read_table(f'{input_dir}/data/evonik_day0vday4_md.tsv', index_col = 0)\n",
+ "metadata = metadata[(metadata['timepoint'] == 'post') & (metadata['chemical'] == 'SLES') & (metadata['individual'] != 'sub10')]\n",
+ "\n",
+ "#load samples\n",
+ "table = load_table(f'{input_dir}/data/day0vday4_raw.biom').to_dataframe().T\n",
+ "table = table.loc[metadata.index.tolist(), syncom_strains]\n",
+ "\n",
+ "# table normalization by column\n",
+ "table = table.loc[:, table.sum(axis=0)>0]\n",
+ "table = table.apply(lambda x: x / x.sum(), axis=1)\n",
+ "\n",
+ "reference_species = 'Corynebacterium afermentans'\n",
+ "\n",
+ "def balance_f(x):\n",
+ " lr = [np.log(x['Cutibacterium acnes']+0.001) - np.log(x[reference_species]+0.001),\n",
+ " np.log(x['Micrococcus luteus']+0.001) - np.log(x[reference_species]+0.001),\n",
+ " np.log(x['Staphylococcus aureus']+0.001) - np.log(x[reference_species]+0.001),\n",
+ " np.log(x['Staphylococcus capitis']+0.001) - np.log(x[reference_species]+0.001),\n",
+ " np.log(x['Staphylococcus epidermidis']+0.001) - np.log(x[reference_species]+0.001),\n",
+ " np.log(x['Staphylococcus hominis']+0.001) - np.log(x[reference_species]+0.001),\n",
+ " np.log(x['Staphylococcus warneri']+0.001) - np.log(x[reference_species]+0.001),\n",
+ " np.log(x['Streptococcus mitis']+0.001) - np.log(x[reference_species]+0.001)\n",
+ " ]\n",
+ " cols = [\n",
+ " 'log(C.acnes/C.afermentans)',\n",
+ " 'log(M.luteus/C.afermentans)',\n",
+ " 'log(S.aureus/C.afermentans)',\n",
+ " 'log(S.capitis/C.afermentans)',\n",
+ " 'log(S.epidermidis/C.afermentans)',\n",
+ " 'log(S.hominis/C.afermentans)',\n",
+ " 'log(S.warneri/C.afermentans)',\n",
+ " 'log(S.mitis/C.afermentans)'\n",
+ " ]\n",
+ " return pd.Series(lr, index=cols)\n",
+ "\n",
+ "balances = table.apply(balance_f, axis = 1)\n",
+ "\n",
+ "f = lambda x: paired_t_test(category='treatment', md=metadata, df=x)[0]\n",
+ "tstats = balances.apply(f, axis=0)\n",
+ "\n",
+ "f = lambda x: paired_t_test(category='treatment', md=metadata, df=x)[1]\n",
+ "tpvals = balances.apply(f, axis=0)\n",
+ "\n",
+ "f = lambda x: paired_t_test(category='treatment', md=metadata, df=x)[2]\n",
+ "tlower = balances.apply(f, axis=0)\n",
+ "\n",
+ "f = lambda x: paired_t_test(category='treatment', md=metadata, df=x)[3]\n",
+ "tupper = balances.apply(f, axis=0)\n",
+ "\n",
+ "stats = pd.DataFrame(\n",
+ " {\n",
+ " '24hControl_vs_24hChemical_tstat': tstats,\n",
+ " '24hControl__vs_24hChemical_pvalue': tpvals,\n",
+ " '24hControl_vs_24hChemical_upperCI': tlower,\n",
+ " '24hControl_vs_24hChemical_lowerCI': tupper\n",
+ " }\n",
+ ")\n",
+ "\n",
+ "data = pd.merge(balances, metadata, left_index=True, right_index=True)\n",
+ "\n",
+ "subdata = data[[\n",
+ " 'log(C.acnes/C.afermentans)',\n",
+ " 'log(M.luteus/C.afermentans)',\n",
+ " 'log(S.aureus/C.afermentans)',\n",
+ " 'log(S.capitis/C.afermentans)',\n",
+ " 'log(S.epidermidis/C.afermentans)',\n",
+ " 'log(S.hominis/C.afermentans)',\n",
+ " 'log(S.warneri/C.afermentans)',\n",
+ " 'log(S.mitis/C.afermentans)',\n",
+ " 'treatment']]\n",
+ "\n",
+ "df = pd.melt(subdata, id_vars = 'treatment', var_name = 'balance')\n",
+ "\n",
+ "# get the palette color name\n",
+ "deep_palette = sns.color_palette(\"Set1\")\n",
+ "\n",
+ "# Extract the first two colors\n",
+ "color1 = deep_palette[0]\n",
+ "color2 = deep_palette[1]\n",
+ "\n",
+ "#, notch=True\n",
+ "sns.set(style=\"white\", font_scale=1.5, font=\"sans-serif\", palette=\"Set1\")\n",
+ "\n",
+ "# get a label for outlier,figure out why there is error bar \n",
+ "# customize which counts as outlier\n",
+ "g = sns.catplot(x='value', y='balance', data=df, hue='treatment', kind='box',\n",
+ " hue_order=['left', 'right'], n_boot = 1000,\n",
+ " height = 10, aspect = 1.8, legend = False, whis = [0,100])\n",
+ "g.set_xlabels('log fold change', fontsize=14)\n",
+ "g.set_ylabels('', fontsize=14)\n",
+ "plt.legend(fontsize = 20, frameon=False)\n",
+ "# plt.legend(loc='lower right', labels=['Control', 'Treatment'])\n",
+ "plt.legend(fontsize=20, frameon=False, loc='lower right')\n",
+ "plt.xticks(fontsize=12)\n",
+ "\n",
+ "ax = plt.gca()\n",
+ "leg = ax.legend(fontsize=20, frameon=False, loc='lower right')\n",
+ "leg.legendHandles[0].set_linewidth(8)\n",
+ "leg.legendHandles[0].set_color(color1)\n",
+ "leg.legendHandles[1].set_linewidth(8)\n",
+ "leg.legendHandles[1].set_color(color2)\n",
+ "leg.get_frame().set_linewidth(0)\n",
+ "#plt.xticks(rotation=60)\n",
+ "\n",
+ "# g.savefig('./out_sles/fig_sles.svg', format='svg')\n",
+ "# stats.to_csv('./out_c_afermentans/stats_day4_sles.csv')\n",
+ "stats"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "ef86565d-5578-45fd-8423-4ad797746748",
+ "metadata": {},
+ "source": [
+ "# RL"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 34,
+ "id": "b7fee009-bdaf-49c7-8eb7-2cc924baa06f",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def paired_t_test(category, md, df):\n",
+ " \"\"\" Performs paired t-test on a dataframe with complementing\n",
+ " sample metadata.\n",
+ "\n",
+ " Parameters\n",
+ " ----------\n",
+ " category : str\n",
+ " Category to test in the dataframe df\n",
+ " md : pd.DataFrame\n",
+ " Sample metadata specifying groupings\n",
+ " df : pd.Dataframe\n",
+ " Counts of microbes across all samples\n",
+ " \n",
+ " Returns\n",
+ " -------\n",
+ " T : float\n",
+ " T-statistic\n",
+ " prob : float\n",
+ " p-value\n",
+ " lower : float\n",
+ " lower 95% CI bound\n",
+ " upper : float\n",
+ " upper 95% CI bound\n",
+ " ddof : int\n",
+ " degrees of freedom\n",
+ " \"\"\"\n",
+ " cats = list(md[category].value_counts().index)\n",
+ " src = md[md[category] == cats[1]].index\n",
+ " print('***')\n",
+ " print(src)\n",
+ " dest = md[md[category] == cats[0]].index\n",
+ " return ttest_rel(df.loc[dest].values, df.loc[src].values)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 35,
+ "id": "43dcf082-3703-40f0-b5af-ca7175ace397",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# load md\n",
+ "metadata = pd.read_table(f'{input_dir}/data/evonik_day0vday4_md.tsv', index_col = 0)\n",
+ "metadata = metadata[(metadata['timepoint'] == 'post') & (metadata['chemical'] == 'RL') & (metadata['individual'] != 'sub4')]\n",
+ "\n",
+ "#load samples\n",
+ "table = load_table(f'{input_dir}/data/day0vday4_raw.biom').to_dataframe().T\n",
+ "table = table.loc[metadata.index.tolist(), syncom_strains]"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 36,
+ "id": "76e7e7d3-a42d-4215-9ce9-579ba02a51c7",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "***\n",
+ "Index(['post_left_sub5', 'post_left_sub6', 'post_left_sub7', 'post_left_sub8',\n",
+ " 'post_left_sub9', 'post_left_sub11', 'post_left_sub3'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['post_left_sub5', 'post_left_sub6', 'post_left_sub7', 'post_left_sub8',\n",
+ " 'post_left_sub9', 'post_left_sub11', 'post_left_sub3'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['post_left_sub5', 'post_left_sub6', 'post_left_sub7', 'post_left_sub8',\n",
+ " 'post_left_sub9', 'post_left_sub11', 'post_left_sub3'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['post_left_sub5', 'post_left_sub6', 'post_left_sub7', 'post_left_sub8',\n",
+ " 'post_left_sub9', 'post_left_sub11', 'post_left_sub3'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['post_left_sub5', 'post_left_sub6', 'post_left_sub7', 'post_left_sub8',\n",
+ " 'post_left_sub9', 'post_left_sub11', 'post_left_sub3'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['post_left_sub5', 'post_left_sub6', 'post_left_sub7', 'post_left_sub8',\n",
+ " 'post_left_sub9', 'post_left_sub11', 'post_left_sub3'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['post_left_sub5', 'post_left_sub6', 'post_left_sub7', 'post_left_sub8',\n",
+ " 'post_left_sub9', 'post_left_sub11', 'post_left_sub3'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['post_left_sub5', 'post_left_sub6', 'post_left_sub7', 'post_left_sub8',\n",
+ " 'post_left_sub9', 'post_left_sub11', 'post_left_sub3'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['post_left_sub5', 'post_left_sub6', 'post_left_sub7', 'post_left_sub8',\n",
+ " 'post_left_sub9', 'post_left_sub11', 'post_left_sub3'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['post_left_sub5', 'post_left_sub6', 'post_left_sub7', 'post_left_sub8',\n",
+ " 'post_left_sub9', 'post_left_sub11', 'post_left_sub3'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['post_left_sub5', 'post_left_sub6', 'post_left_sub7', 'post_left_sub8',\n",
+ " 'post_left_sub9', 'post_left_sub11', 'post_left_sub3'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['post_left_sub5', 'post_left_sub6', 'post_left_sub7', 'post_left_sub8',\n",
+ " 'post_left_sub9', 'post_left_sub11', 'post_left_sub3'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['post_left_sub5', 'post_left_sub6', 'post_left_sub7', 'post_left_sub8',\n",
+ " 'post_left_sub9', 'post_left_sub11', 'post_left_sub3'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['post_left_sub5', 'post_left_sub6', 'post_left_sub7', 'post_left_sub8',\n",
+ " 'post_left_sub9', 'post_left_sub11', 'post_left_sub3'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['post_left_sub5', 'post_left_sub6', 'post_left_sub7', 'post_left_sub8',\n",
+ " 'post_left_sub9', 'post_left_sub11', 'post_left_sub3'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['post_left_sub5', 'post_left_sub6', 'post_left_sub7', 'post_left_sub8',\n",
+ " 'post_left_sub9', 'post_left_sub11', 'post_left_sub3'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['post_left_sub5', 'post_left_sub6', 'post_left_sub7', 'post_left_sub8',\n",
+ " 'post_left_sub9', 'post_left_sub11', 'post_left_sub3'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['post_left_sub5', 'post_left_sub6', 'post_left_sub7', 'post_left_sub8',\n",
+ " 'post_left_sub9', 'post_left_sub11', 'post_left_sub3'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['post_left_sub5', 'post_left_sub6', 'post_left_sub7', 'post_left_sub8',\n",
+ " 'post_left_sub9', 'post_left_sub11', 'post_left_sub3'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['post_left_sub5', 'post_left_sub6', 'post_left_sub7', 'post_left_sub8',\n",
+ " 'post_left_sub9', 'post_left_sub11', 'post_left_sub3'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['post_left_sub5', 'post_left_sub6', 'post_left_sub7', 'post_left_sub8',\n",
+ " 'post_left_sub9', 'post_left_sub11', 'post_left_sub3'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['post_left_sub5', 'post_left_sub6', 'post_left_sub7', 'post_left_sub8',\n",
+ " 'post_left_sub9', 'post_left_sub11', 'post_left_sub3'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['post_left_sub5', 'post_left_sub6', 'post_left_sub7', 'post_left_sub8',\n",
+ " 'post_left_sub9', 'post_left_sub11', 'post_left_sub3'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['post_left_sub5', 'post_left_sub6', 'post_left_sub7', 'post_left_sub8',\n",
+ " 'post_left_sub9', 'post_left_sub11', 'post_left_sub3'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['post_left_sub5', 'post_left_sub6', 'post_left_sub7', 'post_left_sub8',\n",
+ " 'post_left_sub9', 'post_left_sub11', 'post_left_sub3'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['post_left_sub5', 'post_left_sub6', 'post_left_sub7', 'post_left_sub8',\n",
+ " 'post_left_sub9', 'post_left_sub11', 'post_left_sub3'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['post_left_sub5', 'post_left_sub6', 'post_left_sub7', 'post_left_sub8',\n",
+ " 'post_left_sub9', 'post_left_sub11', 'post_left_sub3'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['post_left_sub5', 'post_left_sub6', 'post_left_sub7', 'post_left_sub8',\n",
+ " 'post_left_sub9', 'post_left_sub11', 'post_left_sub3'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['post_left_sub5', 'post_left_sub6', 'post_left_sub7', 'post_left_sub8',\n",
+ " 'post_left_sub9', 'post_left_sub11', 'post_left_sub3'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['post_left_sub5', 'post_left_sub6', 'post_left_sub7', 'post_left_sub8',\n",
+ " 'post_left_sub9', 'post_left_sub11', 'post_left_sub3'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['post_left_sub5', 'post_left_sub6', 'post_left_sub7', 'post_left_sub8',\n",
+ " 'post_left_sub9', 'post_left_sub11', 'post_left_sub3'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['post_left_sub5', 'post_left_sub6', 'post_left_sub7', 'post_left_sub8',\n",
+ " 'post_left_sub9', 'post_left_sub11', 'post_left_sub3'],\n",
+ " dtype='object', name='sampleid')\n"
+ ]
+ },
+ {
+ "data": {
+ "text/html": [
+ "\n",
+ "\n",
+ "
\n",
+ " \n",
+ " \n",
+ " | \n",
+ " 24hControl_vs_24hChemical_tstat | \n",
+ " 24hControl__vs_24hChemical_pvalue | \n",
+ " 24hControl_vs_24hChemical_upperCI | \n",
+ " 24hControl_vs_24hChemical_lowerCI | \n",
+ "
\n",
+ " \n",
+ " \n",
+ " \n",
+ " | log(C.acnes/C.afermentans) | \n",
+ " -0.453956 | \n",
+ " 0.665819 | \n",
+ " -0.912310 | \n",
+ " 0.626776 | \n",
+ "
\n",
+ " \n",
+ " | log(M.luteus/C.afermentans) | \n",
+ " 0.876440 | \n",
+ " 0.414499 | \n",
+ " -0.365185 | \n",
+ " 0.772787 | \n",
+ "
\n",
+ " \n",
+ " | log(S.aureus/C.afermentans) | \n",
+ " -0.273301 | \n",
+ " 0.793786 | \n",
+ " -0.603805 | \n",
+ " 0.482476 | \n",
+ "
\n",
+ " \n",
+ " | log(S.capitis/C.afermentans) | \n",
+ " 0.715330 | \n",
+ " 0.501279 | \n",
+ " -0.464750 | \n",
+ " 0.848735 | \n",
+ "
\n",
+ " \n",
+ " | log(S.epidermidis/C.afermentans) | \n",
+ " 0.364150 | \n",
+ " 0.728235 | \n",
+ " -1.049873 | \n",
+ " 1.416992 | \n",
+ "
\n",
+ " \n",
+ " | log(S.hominis/C.afermentans) | \n",
+ " 0.019290 | \n",
+ " 0.985235 | \n",
+ " -0.766951 | \n",
+ " 0.779140 | \n",
+ "
\n",
+ " \n",
+ " | log(S.warneri/C.afermentans) | \n",
+ " -0.142990 | \n",
+ " 0.890980 | \n",
+ " -0.683197 | \n",
+ " 0.607758 | \n",
+ "
\n",
+ " \n",
+ " | log(S.mitis/C.afermentans) | \n",
+ " 0.258699 | \n",
+ " 0.804514 | \n",
+ " -0.785381 | \n",
+ " 0.971083 | \n",
+ "
\n",
+ " \n",
+ "
\n",
+ "
"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "# table normalization by column\n",
+ "table = table.loc[:, table.sum(axis=0)>0]\n",
+ "table = table.apply(lambda x: x / x.sum(), axis=1)\n",
+ "\n",
+ "reference_species = 'Corynebacterium afermentans'\n",
+ "\n",
+ "def balance_f(x):\n",
+ " lr = [np.log(x['Cutibacterium acnes']+0.001) - np.log(x[reference_species]+0.001),\n",
+ " np.log(x['Micrococcus luteus']+0.001) - np.log(x[reference_species]+0.001),\n",
+ " np.log(x['Staphylococcus aureus']+0.001) - np.log(x[reference_species]+0.001),\n",
+ " np.log(x['Staphylococcus capitis']+0.001) - np.log(x[reference_species]+0.001),\n",
+ " np.log(x['Staphylococcus epidermidis']+0.001) - np.log(x[reference_species]+0.001),\n",
+ " np.log(x['Staphylococcus hominis']+0.001) - np.log(x[reference_species]+0.001),\n",
+ " np.log(x['Staphylococcus warneri']+0.001) - np.log(x[reference_species]+0.001),\n",
+ " np.log(x['Streptococcus mitis']+0.001) - np.log(x[reference_species]+0.001)\n",
+ " ]\n",
+ " cols = [\n",
+ " 'log(C.acnes/C.afermentans)',\n",
+ " 'log(M.luteus/C.afermentans)',\n",
+ " 'log(S.aureus/C.afermentans)',\n",
+ " 'log(S.capitis/C.afermentans)',\n",
+ " 'log(S.epidermidis/C.afermentans)',\n",
+ " 'log(S.hominis/C.afermentans)',\n",
+ " 'log(S.warneri/C.afermentans)',\n",
+ " 'log(S.mitis/C.afermentans)'\n",
+ " ]\n",
+ " return pd.Series(lr, index=cols)\n",
+ "\n",
+ "balances = table.apply(balance_f, axis = 1)\n",
+ "\n",
+ "f = lambda x: paired_t_test(category='treatment', md=metadata, df=x)[0]\n",
+ "tstats = balances.apply(f, axis=0)\n",
+ "\n",
+ "f = lambda x: paired_t_test(category='treatment', md=metadata, df=x)[1]\n",
+ "tpvals = balances.apply(f, axis=0)\n",
+ "\n",
+ "f = lambda x: paired_t_test(category='treatment', md=metadata, df=x)[2]\n",
+ "tlower = balances.apply(f, axis=0)\n",
+ "\n",
+ "f = lambda x: paired_t_test(category='treatment', md=metadata, df=x)[3]\n",
+ "tupper = balances.apply(f, axis=0)\n",
+ "\n",
+ "stats = pd.DataFrame(\n",
+ " {\n",
+ " '24hControl_vs_24hChemical_tstat': tstats,\n",
+ " '24hControl__vs_24hChemical_pvalue': tpvals,\n",
+ " '24hControl_vs_24hChemical_upperCI': tlower,\n",
+ " '24hControl_vs_24hChemical_lowerCI': tupper\n",
+ " }\n",
+ ")\n",
+ "\n",
+ "data = pd.merge(balances, metadata, left_index=True, right_index=True)\n",
+ "\n",
+ "subdata = data[[\n",
+ " 'log(C.acnes/C.afermentans)',\n",
+ " 'log(M.luteus/C.afermentans)',\n",
+ " 'log(S.aureus/C.afermentans)',\n",
+ " 'log(S.capitis/C.afermentans)',\n",
+ " 'log(S.epidermidis/C.afermentans)',\n",
+ " 'log(S.hominis/C.afermentans)',\n",
+ " 'log(S.warneri/C.afermentans)',\n",
+ " 'log(S.mitis/C.afermentans)',\n",
+ " 'treatment']]\n",
+ "\n",
+ "df = pd.melt(subdata, id_vars = 'treatment', var_name = 'balance')\n",
+ "\n",
+ "# get the palette color name\n",
+ "deep_palette = sns.color_palette(\"Set1\")\n",
+ "\n",
+ "# Extract the first two colors\n",
+ "color1 = deep_palette[0]\n",
+ "color2 = deep_palette[1]\n",
+ "\n",
+ "#, notch=True\n",
+ "sns.set(style=\"white\", font_scale=1.5, font=\"sans-serif\", palette=\"Set1\")\n",
+ "\n",
+ "# get a label for outlier,figure out why there is error bar \n",
+ "# customize which counts as outlier\n",
+ "g = sns.catplot(x='value', y='balance', data=df, hue='treatment', kind='box',\n",
+ " hue_order=['left', 'right'], n_boot = 1000,\n",
+ " height = 10, aspect = 1.8, legend = False, whis = [0,100])\n",
+ "g.set_xlabels('log fold change', fontsize=14)\n",
+ "g.set_ylabels('', fontsize=14)\n",
+ "plt.legend(fontsize = 20, frameon=False)\n",
+ "# plt.legend(loc='lower right', labels=['Control', 'Treatment'])\n",
+ "plt.legend(fontsize=20, frameon=False, loc='lower right')\n",
+ "plt.xticks(fontsize=12)\n",
+ "\n",
+ "ax = plt.gca()\n",
+ "leg = ax.legend(fontsize=20, frameon=False, loc='lower right')\n",
+ "leg.legendHandles[0].set_linewidth(8)\n",
+ "leg.legendHandles[0].set_color(color1)\n",
+ "leg.legendHandles[1].set_linewidth(8)\n",
+ "leg.legendHandles[1].set_color(color2)\n",
+ "leg.get_frame().set_linewidth(0)\n",
+ "#plt.xticks(rotation=60)\n",
+ "\n",
+ "# g.savefig('./out_sles/fig_sles.svg', format='svg')\n",
+ "# stats.to_csv('./out_c_afermentans/stats_day4_rl.csv')\n",
+ "stats"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "9e1cacb4-26dd-4c18-b5a8-dae7a8d75bb3",
+ "metadata": {},
+ "outputs": [],
+ "source": []
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "qiime2-2020.6-system",
+ "language": "python",
+ "name": "qiime2-2020.6-system"
+ },
+ "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.10"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/figure_generation_clinical_day7_afermentans_reordered.ipynb b/figure_generation_clinical_day7_afermentans_reordered.ipynb
new file mode 100644
index 0000000..495852a
--- /dev/null
+++ b/figure_generation_clinical_day7_afermentans_reordered.ipynb
@@ -0,0 +1,1292 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "id": "3f300489-27ee-4dd4-bfed-bae15a09a71f",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# ----------------------------------------------------------------------------\n",
+ "# SkinCom Analysis - Perform t-test & boxplots for log-ratios\n",
+ "# Author: Sherlyn Weng\n",
+ "# Version: 1.0.0\n",
+ "# Mmaintainer: Sherlyn Weng\n",
+ "# Email: y1weng@ucsd.edu\n",
+ "# This code is adapted from https://github.com/knightlab-analyses/reference-frames\n",
+ "# ----------------------------------------------------------------------------"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "id": "201632e3-64c2-4f20-88ee-9166b9760e4e",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import numpy as np\n",
+ "import pandas as pd\n",
+ "import seaborn as sns\n",
+ "import matplotlib.pyplot as plt\n",
+ "import matplotlib as mpl\n",
+ "from biom import load_table\n",
+ "from scipy.stats import ttest_rel, wilcoxon, ttest_ind, pearsonr, spearmanr\n",
+ "from skbio.stats.composition import multiplicative_replacement\n",
+ "import matplotlib.gridspec as gridspec\n",
+ "\n",
+ "from collections import namedtuple\n",
+ "\n",
+ "%matplotlib inline\n",
+ "\n",
+ "import warnings\n",
+ "warnings.filterwarnings('ignore')\n",
+ "\n",
+ "\"\"\"\n",
+ "Here, we have implemented an extended version of scipy's paired t-test function.\n",
+ "see https://github.com/scipy/scipy/blob/master/scipy/stats/mstats_basic.py#L1085\n",
+ "\n",
+ "This should really be pushed to scipy -- any volunteers?\n",
+ "\"\"\"\n",
+ "import numpy as np\n",
+ "from scipy.stats import t\n",
+ "from scipy.special import betainc\n",
+ "from collections import namedtuple\n",
+ "\n",
+ "\n",
+ "Ttest_relResult = namedtuple('Ttest_relResult', \n",
+ " ('statistic', 'pvalue', 'lower_CI', 'upper_CI', 'ddof'))\n",
+ "\n",
+ "\n",
+ "def ttest_rel(x, y):\n",
+ " \n",
+ " D = x - y\n",
+ " n = len(D)\n",
+ " \n",
+ " sd = np.std(D, ddof=1) \n",
+ " se = sd / np.sqrt(n)\n",
+ " mu = np.mean(D)\n",
+ " T = mu / se\n",
+ " ddof = n-1\n",
+ " ci = se * t.ppf( 0.05 / 2, ddof) \n",
+ " lower = mu - np.abs(ci)\n",
+ " upper = mu + np.abs(ci)\n",
+ " \n",
+ " prob = betainc(0.5*ddof, 0.5, ddof/(ddof + T*T))\n",
+ " return Ttest_relResult(T, prob, lower, upper, int(ddof))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 6,
+ "id": "2e1793ae-1a5d-4bc1-a63a-2beb2ea0dd54",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def paired_t_test(category, md, df):\n",
+ " \"\"\" Performs paired t-test on a dataframe with complementing\n",
+ " sample metadata.\n",
+ "\n",
+ " Parameters\n",
+ " ----------\n",
+ " category : str\n",
+ " Category to test in the dataframe df\n",
+ " md : pd.DataFrame\n",
+ " Sample metadata specifying groupings\n",
+ " df : pd.Dataframe\n",
+ " Counts of microbes across all samples\n",
+ " \n",
+ " Returns\n",
+ " -------\n",
+ " T : float\n",
+ " T-statistic\n",
+ " prob : float\n",
+ " p-value\n",
+ " lower : float\n",
+ " lower 95% CI bound\n",
+ " upper : float\n",
+ " upper 95% CI bound\n",
+ " ddof : int\n",
+ " degrees of freedom\n",
+ " \"\"\"\n",
+ " cats = list(md[category].value_counts().index)\n",
+ " src = md[md[category] == cats[0]].index\n",
+ " print('hio')\n",
+ " print(src)\n",
+ " dest = md[md[category] == cats[1]].index\n",
+ " return ttest_rel(df.loc[dest].values, df.loc[src].values)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "5242c6a2-b6b2-497e-aa34-225f3a9c1cab",
+ "metadata": {},
+ "source": [
+ "## Input Directory"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 7,
+ "id": "0f656bb7-fbf6-4b0e-909b-d94835646df4",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "input_dir = \"/home/y1weng/02_evonik_batch2/songbird/evonik\""
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "f9b78015-1a9a-4d8e-97d0-e7ff85a1e82c",
+ "metadata": {},
+ "source": [
+ "## Variables"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "16391819-2ad3-4cd3-b051-0abb8d737f98",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# # load metadata\n",
+ "# metadata = pd.read_table(f'{input_dir}/data/evonik_day0vday4_md.tsv', index_col = 0)\n",
+ "# filtered_metadata = metadata[(metadata['time_treatment'].isin(['post_control', 'post_RL', 'post_SLES'])) &\n",
+ "# (metadata['individual'].isin(['sub1', 'sub2','sub10', 'sub11', 'sub12', 'sub13', 'sub5', 'sub6', 'sub8', 'sub9'])) \n",
+ "# ].sort_values('individual')\n",
+ "# filtered_metadata_rl = metadata[(metadata['time_treatment'].isin(['post_control', 'post_RL', 'post_SLES'])) &\n",
+ "# (metadata['individual'].isin(['sub1', 'sub2', 'sub10', 'sub11', 'sub12', 'sub13', 'sub5', 'sub6', 'sub8', 'sub9'])) & \n",
+ "# (metadata['individual'].isin(['sub3', 'sub4', 'sub5', 'sub6', 'sub7', 'sub8', 'sub9', 'sub11'])) \n",
+ "# ].sort_values('individual')\n",
+ "# filtered_metadata_sles = metadata[(metadata['time_treatment'].isin(['post_control', 'post_RL', 'post_SLES'])) &\n",
+ "# (metadata['individual'].isin(['sub1', 'sub2','sub10', 'sub11', 'sub12', 'sub13', 'sub5', 'sub6', 'sub8', 'sub9'])) &\n",
+ "# (metadata['individual'].isin(['sub1', 'sub2', 'sub12', 'sub13'])) \n",
+ "# ].sort_values('individual')\n",
+ "# #metadata = filtered_metadata_sles\n",
+ "# #metadata = filtered_metadata"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 8,
+ "id": "af3b610d-80e5-4525-b1bf-86d36a88e9cb",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "syncom_strains = [\"Corynebacterium afermentans\",\n",
+ " \"Cutibacterium acnes\",\n",
+ " \"Micrococcus luteus\",\n",
+ " \"Staphylococcus aureus\",\n",
+ " \"Staphylococcus capitis\",\n",
+ " \"Staphylococcus epidermidis\",\n",
+ " \"Staphylococcus hominis\",\n",
+ " \"Staphylococcus warneri\",\n",
+ " \"Streptococcus mitis\"]"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 49,
+ "id": "8e5579b5-465f-4558-bc6f-4a0698f8dcf7",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# # modify metadata\n",
+ "# metadata = pd.read_table(f'{input_dir}/evonik_songbird_metadata.txt', index_col = 0)\n",
+ "# metadata['chemical'] = np.where(metadata['Subject'].isin(['sub1', 'sub2', 'sub10', 'sub12', 'sub13']), 'SLES', 'RL')\n",
+ "\n",
+ "# def determine_time_treatment(row):\n",
+ "# if row['Time_Point'] == 'pre':\n",
+ "# return 'pre_control'\n",
+ "# elif row['Time_Point'] == 'post' and row['Control'] == 'left':\n",
+ "# return 'post_control'\n",
+ "# elif row['Time_Point'] == 'post' and row['Control'] == 'right' and row['chemical'] == 'SLES':\n",
+ "# return 'post_SLES'\n",
+ "# elif row['Time_Point'] == 'post' and row['Control'] == 'right' and row['chemical'] == 'RL':\n",
+ "# return 'post_RL'\n",
+ "# else:\n",
+ "# return 'unknown'\n",
+ "\n",
+ "# metadata['time_treatment'] = metadata.apply(determine_time_treatment, axis = 1)\n",
+ "\n",
+ "# metadata = metadata.sort_values('time_treatment')\n",
+ "\n",
+ "# testing_list = ['Train'] * 3 + ['Test'] * 2 + ['Train'] * 3 + ['Test'] * 2 + ['Train'] * 9 + ['Test'] * 3 + ['Train'] * 13 + ['Test'] * 4\n",
+ "\n",
+ "# metadata['testing'] = testing_list\n",
+ "\n",
+ "# metadata\n",
+ "\n",
+ "# metadata.to_csv('/home/y1weng/02_evonik_batch2/songbird/evonik/day7_metadata.tsv', sep='\\t', index=True, header=True)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 47,
+ "id": "2406bf00-ebc3-4453-bcda-5170212f866b",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# # create a feature md\n",
+ "# df = pd.read_csv('/home/y1weng/02_evonik_batch2/songbird/evonik/evonik_all_taxonomy_reformatted.tsv', sep='\\t')\n",
+ "# df['sampleid'] = df.featureid\n",
+ "# df['full_name'] = df.featureid\n",
+ "\n",
+ "# df = df.loc[:, ['sampleid', 'full_name']]\n",
+ "# df.to_csv('/home/y1weng/02_evonik_batch2/songbird/evonik/day7_featuremd.tsv', sep='\\t', index=False, header=True)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "1164688b-1357-4d73-bedb-66dc5b71386a",
+ "metadata": {
+ "tags": []
+ },
+ "source": [
+ "## SLES"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 14,
+ "id": "d7fc69b7-01f2-449a-ad95-f7cb6dd60cda",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "hio\n",
+ "Index(['post_left_sub1_S107_WoL_index_filtered',\n",
+ " 'post_left_sub10_S15_WoL_index_filtered',\n",
+ " 'post_left_sub12_S17_WoL_index_filtered',\n",
+ " 'post_left_sub13_S111_WoL_index_filtered'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hio\n",
+ "Index(['post_left_sub1_S107_WoL_index_filtered',\n",
+ " 'post_left_sub10_S15_WoL_index_filtered',\n",
+ " 'post_left_sub12_S17_WoL_index_filtered',\n",
+ " 'post_left_sub13_S111_WoL_index_filtered'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hio\n",
+ "Index(['post_left_sub1_S107_WoL_index_filtered',\n",
+ " 'post_left_sub10_S15_WoL_index_filtered',\n",
+ " 'post_left_sub12_S17_WoL_index_filtered',\n",
+ " 'post_left_sub13_S111_WoL_index_filtered'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hio\n",
+ "Index(['post_left_sub1_S107_WoL_index_filtered',\n",
+ " 'post_left_sub10_S15_WoL_index_filtered',\n",
+ " 'post_left_sub12_S17_WoL_index_filtered',\n",
+ " 'post_left_sub13_S111_WoL_index_filtered'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hio\n",
+ "Index(['post_left_sub1_S107_WoL_index_filtered',\n",
+ " 'post_left_sub10_S15_WoL_index_filtered',\n",
+ " 'post_left_sub12_S17_WoL_index_filtered',\n",
+ " 'post_left_sub13_S111_WoL_index_filtered'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hio\n",
+ "Index(['post_left_sub1_S107_WoL_index_filtered',\n",
+ " 'post_left_sub10_S15_WoL_index_filtered',\n",
+ " 'post_left_sub12_S17_WoL_index_filtered',\n",
+ " 'post_left_sub13_S111_WoL_index_filtered'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hio\n",
+ "Index(['post_left_sub1_S107_WoL_index_filtered',\n",
+ " 'post_left_sub10_S15_WoL_index_filtered',\n",
+ " 'post_left_sub12_S17_WoL_index_filtered',\n",
+ " 'post_left_sub13_S111_WoL_index_filtered'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hio\n",
+ "Index(['post_left_sub1_S107_WoL_index_filtered',\n",
+ " 'post_left_sub10_S15_WoL_index_filtered',\n",
+ " 'post_left_sub12_S17_WoL_index_filtered',\n",
+ " 'post_left_sub13_S111_WoL_index_filtered'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hio\n",
+ "Index(['post_left_sub1_S107_WoL_index_filtered',\n",
+ " 'post_left_sub10_S15_WoL_index_filtered',\n",
+ " 'post_left_sub12_S17_WoL_index_filtered',\n",
+ " 'post_left_sub13_S111_WoL_index_filtered'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hio\n",
+ "Index(['post_left_sub1_S107_WoL_index_filtered',\n",
+ " 'post_left_sub10_S15_WoL_index_filtered',\n",
+ " 'post_left_sub12_S17_WoL_index_filtered',\n",
+ " 'post_left_sub13_S111_WoL_index_filtered'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hio\n",
+ "Index(['post_left_sub1_S107_WoL_index_filtered',\n",
+ " 'post_left_sub10_S15_WoL_index_filtered',\n",
+ " 'post_left_sub12_S17_WoL_index_filtered',\n",
+ " 'post_left_sub13_S111_WoL_index_filtered'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hio\n",
+ "Index(['post_left_sub1_S107_WoL_index_filtered',\n",
+ " 'post_left_sub10_S15_WoL_index_filtered',\n",
+ " 'post_left_sub12_S17_WoL_index_filtered',\n",
+ " 'post_left_sub13_S111_WoL_index_filtered'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hio\n",
+ "Index(['post_left_sub1_S107_WoL_index_filtered',\n",
+ " 'post_left_sub10_S15_WoL_index_filtered',\n",
+ " 'post_left_sub12_S17_WoL_index_filtered',\n",
+ " 'post_left_sub13_S111_WoL_index_filtered'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hio\n",
+ "Index(['post_left_sub1_S107_WoL_index_filtered',\n",
+ " 'post_left_sub10_S15_WoL_index_filtered',\n",
+ " 'post_left_sub12_S17_WoL_index_filtered',\n",
+ " 'post_left_sub13_S111_WoL_index_filtered'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hio\n",
+ "Index(['post_left_sub1_S107_WoL_index_filtered',\n",
+ " 'post_left_sub10_S15_WoL_index_filtered',\n",
+ " 'post_left_sub12_S17_WoL_index_filtered',\n",
+ " 'post_left_sub13_S111_WoL_index_filtered'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hio\n",
+ "Index(['post_left_sub1_S107_WoL_index_filtered',\n",
+ " 'post_left_sub10_S15_WoL_index_filtered',\n",
+ " 'post_left_sub12_S17_WoL_index_filtered',\n",
+ " 'post_left_sub13_S111_WoL_index_filtered'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hio\n",
+ "Index(['post_left_sub1_S107_WoL_index_filtered',\n",
+ " 'post_left_sub10_S15_WoL_index_filtered',\n",
+ " 'post_left_sub12_S17_WoL_index_filtered',\n",
+ " 'post_left_sub13_S111_WoL_index_filtered'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hio\n",
+ "Index(['post_left_sub1_S107_WoL_index_filtered',\n",
+ " 'post_left_sub10_S15_WoL_index_filtered',\n",
+ " 'post_left_sub12_S17_WoL_index_filtered',\n",
+ " 'post_left_sub13_S111_WoL_index_filtered'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hio\n",
+ "Index(['post_left_sub1_S107_WoL_index_filtered',\n",
+ " 'post_left_sub10_S15_WoL_index_filtered',\n",
+ " 'post_left_sub12_S17_WoL_index_filtered',\n",
+ " 'post_left_sub13_S111_WoL_index_filtered'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hio\n",
+ "Index(['post_left_sub1_S107_WoL_index_filtered',\n",
+ " 'post_left_sub10_S15_WoL_index_filtered',\n",
+ " 'post_left_sub12_S17_WoL_index_filtered',\n",
+ " 'post_left_sub13_S111_WoL_index_filtered'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hio\n",
+ "Index(['post_left_sub1_S107_WoL_index_filtered',\n",
+ " 'post_left_sub10_S15_WoL_index_filtered',\n",
+ " 'post_left_sub12_S17_WoL_index_filtered',\n",
+ " 'post_left_sub13_S111_WoL_index_filtered'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hio\n",
+ "Index(['post_left_sub1_S107_WoL_index_filtered',\n",
+ " 'post_left_sub10_S15_WoL_index_filtered',\n",
+ " 'post_left_sub12_S17_WoL_index_filtered',\n",
+ " 'post_left_sub13_S111_WoL_index_filtered'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hio\n",
+ "Index(['post_left_sub1_S107_WoL_index_filtered',\n",
+ " 'post_left_sub10_S15_WoL_index_filtered',\n",
+ " 'post_left_sub12_S17_WoL_index_filtered',\n",
+ " 'post_left_sub13_S111_WoL_index_filtered'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hio\n",
+ "Index(['post_left_sub1_S107_WoL_index_filtered',\n",
+ " 'post_left_sub10_S15_WoL_index_filtered',\n",
+ " 'post_left_sub12_S17_WoL_index_filtered',\n",
+ " 'post_left_sub13_S111_WoL_index_filtered'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hio\n",
+ "Index(['post_left_sub1_S107_WoL_index_filtered',\n",
+ " 'post_left_sub10_S15_WoL_index_filtered',\n",
+ " 'post_left_sub12_S17_WoL_index_filtered',\n",
+ " 'post_left_sub13_S111_WoL_index_filtered'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hio\n",
+ "Index(['post_left_sub1_S107_WoL_index_filtered',\n",
+ " 'post_left_sub10_S15_WoL_index_filtered',\n",
+ " 'post_left_sub12_S17_WoL_index_filtered',\n",
+ " 'post_left_sub13_S111_WoL_index_filtered'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hio\n",
+ "Index(['post_left_sub1_S107_WoL_index_filtered',\n",
+ " 'post_left_sub10_S15_WoL_index_filtered',\n",
+ " 'post_left_sub12_S17_WoL_index_filtered',\n",
+ " 'post_left_sub13_S111_WoL_index_filtered'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hio\n",
+ "Index(['post_left_sub1_S107_WoL_index_filtered',\n",
+ " 'post_left_sub10_S15_WoL_index_filtered',\n",
+ " 'post_left_sub12_S17_WoL_index_filtered',\n",
+ " 'post_left_sub13_S111_WoL_index_filtered'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hio\n",
+ "Index(['post_left_sub1_S107_WoL_index_filtered',\n",
+ " 'post_left_sub10_S15_WoL_index_filtered',\n",
+ " 'post_left_sub12_S17_WoL_index_filtered',\n",
+ " 'post_left_sub13_S111_WoL_index_filtered'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hio\n",
+ "Index(['post_left_sub1_S107_WoL_index_filtered',\n",
+ " 'post_left_sub10_S15_WoL_index_filtered',\n",
+ " 'post_left_sub12_S17_WoL_index_filtered',\n",
+ " 'post_left_sub13_S111_WoL_index_filtered'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hio\n",
+ "Index(['post_left_sub1_S107_WoL_index_filtered',\n",
+ " 'post_left_sub10_S15_WoL_index_filtered',\n",
+ " 'post_left_sub12_S17_WoL_index_filtered',\n",
+ " 'post_left_sub13_S111_WoL_index_filtered'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hio\n",
+ "Index(['post_left_sub1_S107_WoL_index_filtered',\n",
+ " 'post_left_sub10_S15_WoL_index_filtered',\n",
+ " 'post_left_sub12_S17_WoL_index_filtered',\n",
+ " 'post_left_sub13_S111_WoL_index_filtered'],\n",
+ " dtype='object', name='sampleid')\n"
+ ]
+ },
+ {
+ "data": {
+ "text/html": [
+ "\n",
+ "\n",
+ "
\n",
+ " \n",
+ " \n",
+ " | \n",
+ " 24hControl_vs_24hChemical_tstat | \n",
+ " 24hControl__vs_24hChemical_pvalue | \n",
+ " 24hControl_vs_24hChemical_upperCI | \n",
+ " 24hControl_vs_24hChemical_lowerCI | \n",
+ "
\n",
+ " \n",
+ " \n",
+ " \n",
+ " | log(C.acnes/C.afermentans) | \n",
+ " -1.869453 | \n",
+ " 0.158339 | \n",
+ " -0.233951 | \n",
+ " 0.060804 | \n",
+ "
\n",
+ " \n",
+ " | log(M.luteus/C.afermentans) | \n",
+ " 0.888782 | \n",
+ " 0.439591 | \n",
+ " -0.548265 | \n",
+ " 0.973165 | \n",
+ "
\n",
+ " \n",
+ " | log(S.aureus/C.afermentans) | \n",
+ " 1.494316 | \n",
+ " 0.231953 | \n",
+ " -0.177493 | \n",
+ " 0.491722 | \n",
+ "
\n",
+ " \n",
+ " | log(S.capitis/C.afermentans) | \n",
+ " 1.762337 | \n",
+ " 0.176220 | \n",
+ " -0.475287 | \n",
+ " 1.654939 | \n",
+ "
\n",
+ " \n",
+ " | log(S.epidermidis/C.afermentans) | \n",
+ " 3.206298 | \n",
+ " 0.049095 | \n",
+ " 0.007707 | \n",
+ " 2.064444 | \n",
+ "
\n",
+ " \n",
+ " | log(S.hominis/C.afermentans) | \n",
+ " 2.089460 | \n",
+ " 0.127843 | \n",
+ " -0.170947 | \n",
+ " 0.824543 | \n",
+ "
\n",
+ " \n",
+ " | log(S.warneri/C.afermentans) | \n",
+ " 1.337519 | \n",
+ " 0.273415 | \n",
+ " -0.592327 | \n",
+ " 1.451166 | \n",
+ "
\n",
+ " \n",
+ " | log(S.mitis/C.afermentans) | \n",
+ " -0.535378 | \n",
+ " 0.629533 | \n",
+ " -0.216667 | \n",
+ " 0.154266 | \n",
+ "
\n",
+ " \n",
+ "
\n",
+ "
"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "# load md\n",
+ "metadata = pd.read_table(f'{input_dir}/evonik_songbird_metadata.txt', index_col = 0)\n",
+ "metadata['chemical'] = np.where(metadata['Subject'].isin(['sub1', 'sub2', 'sub10', 'sub12', 'sub13']), 'SLES', 'RL')\n",
+ "metadata = metadata.sort_values('Subject')\n",
+ "metadata = metadata[(metadata['Time_Point'] == 'post') & (metadata['chemical'] == 'SLES') & (metadata['Subject'] != 'sub2')]\n",
+ "\n",
+ "#load samples\n",
+ "table = load_table(f'{input_dir}/evonik_all_taxonomy_reformatted.biom').to_dataframe().T\n",
+ "table = table.loc[metadata.index.tolist(), syncom_strains]\n",
+ "\n",
+ "# table normalization by column\n",
+ "table = table.loc[:, table.sum(axis=0)>0]\n",
+ "table = table.apply(lambda x: x / x.sum(), axis=1)\n",
+ "\n",
+ "reference_species = 'Corynebacterium afermentans'\n",
+ "\n",
+ "def balance_f(x):\n",
+ " lr = [np.log(x['Cutibacterium acnes']+0.001) - np.log(x[reference_species]+0.001),\n",
+ " np.log(x['Micrococcus luteus']+0.001) - np.log(x[reference_species]+0.001),\n",
+ " np.log(x['Staphylococcus aureus']+0.001) - np.log(x[reference_species]+0.001),\n",
+ " np.log(x['Staphylococcus capitis']+0.001) - np.log(x[reference_species]+0.001),\n",
+ " np.log(x['Staphylococcus epidermidis']+0.001) - np.log(x[reference_species]+0.001),\n",
+ " np.log(x['Staphylococcus hominis']+0.001) - np.log(x[reference_species]+0.001),\n",
+ " np.log(x['Staphylococcus warneri']+0.001) - np.log(x[reference_species]+0.001),\n",
+ " np.log(x['Streptococcus mitis']+0.001) - np.log(x[reference_species]+0.001)\n",
+ " ]\n",
+ " cols = [\n",
+ " 'log(C.acnes/C.afermentans)',\n",
+ " 'log(M.luteus/C.afermentans)',\n",
+ " 'log(S.aureus/C.afermentans)',\n",
+ " 'log(S.capitis/C.afermentans)',\n",
+ " 'log(S.epidermidis/C.afermentans)',\n",
+ " 'log(S.hominis/C.afermentans)',\n",
+ " 'log(S.warneri/C.afermentans)',\n",
+ " 'log(S.mitis/C.afermentans)'\n",
+ " ]\n",
+ " return pd.Series(lr, index=cols)\n",
+ "\n",
+ "balances = table.apply(balance_f, axis = 1)\n",
+ "\n",
+ "f = lambda x: paired_t_test(category='Control', md=metadata, df=x)[0]\n",
+ "tstats = balances.apply(f, axis=0)\n",
+ "\n",
+ "f = lambda x: paired_t_test(category='Control', md=metadata, df=x)[1]\n",
+ "tpvals = balances.apply(f, axis=0)\n",
+ "\n",
+ "f = lambda x: paired_t_test(category='Control', md=metadata, df=x)[2]\n",
+ "tlower = balances.apply(f, axis=0)\n",
+ "\n",
+ "f = lambda x: paired_t_test(category='Control', md=metadata, df=x)[3]\n",
+ "tupper = balances.apply(f, axis=0)\n",
+ "\n",
+ "stats = pd.DataFrame(\n",
+ " {\n",
+ " '24hControl_vs_24hChemical_tstat': tstats,\n",
+ " '24hControl__vs_24hChemical_pvalue': tpvals,\n",
+ " '24hControl_vs_24hChemical_upperCI': tlower,\n",
+ " '24hControl_vs_24hChemical_lowerCI': tupper\n",
+ " }\n",
+ ")\n",
+ "\n",
+ "data = pd.merge(balances, metadata, left_index=True, right_index=True)\n",
+ "\n",
+ "subdata = data[[\n",
+ " 'log(C.acnes/C.afermentans)',\n",
+ " 'log(M.luteus/C.afermentans)',\n",
+ " 'log(S.aureus/C.afermentans)',\n",
+ " 'log(S.capitis/C.afermentans)',\n",
+ " 'log(S.epidermidis/C.afermentans)',\n",
+ " 'log(S.hominis/C.afermentans)',\n",
+ " 'log(S.warneri/C.afermentans)',\n",
+ " 'log(S.mitis/C.afermentans)',\n",
+ " 'Control']]\n",
+ "\n",
+ "df = pd.melt(subdata, id_vars = 'Control', var_name = 'balance')\n",
+ "\n",
+ "# get the palette color name\n",
+ "deep_palette = sns.color_palette(\"Set1\")\n",
+ "\n",
+ "# Extract the first two colors\n",
+ "color1 = deep_palette[0]\n",
+ "color2 = deep_palette[1]\n",
+ "\n",
+ "#, notch=True\n",
+ "sns.set(style=\"white\", font_scale=1.5, font=\"sans-serif\", palette=\"Set1\")\n",
+ "\n",
+ "# get a label for outlier,figure out why there is error bar \n",
+ "# customize which counts as outlier\n",
+ "g = sns.catplot(x='value', y='balance', data=df, hue='Control', kind='box',\n",
+ " hue_order=['left', 'right'], n_boot = 1000,\n",
+ " height = 10, aspect = 1.8, legend = False, whis = [0,100])\n",
+ "g.set_xlabels('log fold change', fontsize=14)\n",
+ "g.set_ylabels('', fontsize=14)\n",
+ "plt.legend(fontsize = 20, frameon=False)\n",
+ "# plt.legend(loc='lower right', labels=['Control', 'Treatment'])\n",
+ "plt.legend(fontsize=20, frameon=False, loc='lower right')\n",
+ "plt.xticks(fontsize=12)\n",
+ "\n",
+ "ax = plt.gca()\n",
+ "leg = ax.legend(fontsize=20, frameon=False, loc='lower right')\n",
+ "leg.legendHandles[0].set_linewidth(8)\n",
+ "leg.legendHandles[0].set_color(color1)\n",
+ "leg.legendHandles[1].set_linewidth(8)\n",
+ "leg.legendHandles[1].set_color(color2)\n",
+ "leg.get_frame().set_linewidth(0)\n",
+ "#plt.xticks(rotation=60)\n",
+ "\n",
+ "# g.savefig('./out_sles/fig_sles.svg', format='svg')\n",
+ "# stats.to_csv('./out_c_afermentans/stats_day7_sles.csv')\n",
+ "stats"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "ef86565d-5578-45fd-8423-4ad797746748",
+ "metadata": {},
+ "source": [
+ "# RL"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 18,
+ "id": "b2fd95e0-6312-4f4c-929b-e0aa03a29669",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def paired_t_test(category, md, df):\n",
+ " \"\"\" Performs paired t-test on a dataframe with complementing\n",
+ " sample metadata.\n",
+ "\n",
+ " Parameters\n",
+ " ----------\n",
+ " category : str\n",
+ " Category to test in the dataframe df\n",
+ " md : pd.DataFrame\n",
+ " Sample metadata specifying groupings\n",
+ " df : pd.Dataframe\n",
+ " Counts of microbes across all samples\n",
+ " \n",
+ " Returns\n",
+ " -------\n",
+ " T : float\n",
+ " T-statistic\n",
+ " prob : float\n",
+ " p-value\n",
+ " lower : float\n",
+ " lower 95% CI bound\n",
+ " upper : float\n",
+ " upper 95% CI bound\n",
+ " ddof : int\n",
+ " degrees of freedom\n",
+ " \"\"\"\n",
+ " cats = list(md[category].value_counts().index)\n",
+ " src = md[md[category] == cats[1]].index\n",
+ " print('***')\n",
+ " print(src)\n",
+ " dest = md[md[category] == cats[0]].index\n",
+ " return ttest_rel(df.loc[dest].values, df.loc[src].values)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 19,
+ "id": "aa36f30c-0a35-4688-922d-4969f8db38b3",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# load md\n",
+ "metadata = pd.read_table(f'{input_dir}/evonik_songbird_metadata.txt', index_col = 0)\n",
+ "metadata['chemical'] = np.where(metadata['Subject'].isin(['sub1', 'sub2', 'sub10', 'sub12', 'sub13']), 'SLES', 'RL')\n",
+ "metadata = metadata.sort_values('Subject')\n",
+ "metadata = metadata[(metadata['Time_Point'] == 'post') & (metadata['chemical'] == 'RL') & (metadata['Subject'] != 'sub3') & (metadata['Subject'] != 'sub4') & (metadata['Subject'] != 'sub7')]\n",
+ "\n",
+ "#load samples\n",
+ "table = load_table(f'{input_dir}/evonik_all_taxonomy_reformatted.biom').to_dataframe().T\n",
+ "table = table.loc[metadata.index.tolist(), syncom_strains]"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 22,
+ "id": "6d035310-3444-459b-9774-5e78a4242ec1",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "***\n",
+ "Index(['post_left_sub11_S16_WoL_index_filtered',\n",
+ " 'post_left_sub5_S12_WoL_index_filtered',\n",
+ " 'post_left_sub6_S13_WoL_index_filtered',\n",
+ " 'post_left_sub8_S14_WoL_index_filtered',\n",
+ " 'post_left_sub9_S110_WoL_index_filtered'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['post_left_sub11_S16_WoL_index_filtered',\n",
+ " 'post_left_sub5_S12_WoL_index_filtered',\n",
+ " 'post_left_sub6_S13_WoL_index_filtered',\n",
+ " 'post_left_sub8_S14_WoL_index_filtered',\n",
+ " 'post_left_sub9_S110_WoL_index_filtered'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['post_left_sub11_S16_WoL_index_filtered',\n",
+ " 'post_left_sub5_S12_WoL_index_filtered',\n",
+ " 'post_left_sub6_S13_WoL_index_filtered',\n",
+ " 'post_left_sub8_S14_WoL_index_filtered',\n",
+ " 'post_left_sub9_S110_WoL_index_filtered'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['post_left_sub11_S16_WoL_index_filtered',\n",
+ " 'post_left_sub5_S12_WoL_index_filtered',\n",
+ " 'post_left_sub6_S13_WoL_index_filtered',\n",
+ " 'post_left_sub8_S14_WoL_index_filtered',\n",
+ " 'post_left_sub9_S110_WoL_index_filtered'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['post_left_sub11_S16_WoL_index_filtered',\n",
+ " 'post_left_sub5_S12_WoL_index_filtered',\n",
+ " 'post_left_sub6_S13_WoL_index_filtered',\n",
+ " 'post_left_sub8_S14_WoL_index_filtered',\n",
+ " 'post_left_sub9_S110_WoL_index_filtered'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['post_left_sub11_S16_WoL_index_filtered',\n",
+ " 'post_left_sub5_S12_WoL_index_filtered',\n",
+ " 'post_left_sub6_S13_WoL_index_filtered',\n",
+ " 'post_left_sub8_S14_WoL_index_filtered',\n",
+ " 'post_left_sub9_S110_WoL_index_filtered'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['post_left_sub11_S16_WoL_index_filtered',\n",
+ " 'post_left_sub5_S12_WoL_index_filtered',\n",
+ " 'post_left_sub6_S13_WoL_index_filtered',\n",
+ " 'post_left_sub8_S14_WoL_index_filtered',\n",
+ " 'post_left_sub9_S110_WoL_index_filtered'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['post_left_sub11_S16_WoL_index_filtered',\n",
+ " 'post_left_sub5_S12_WoL_index_filtered',\n",
+ " 'post_left_sub6_S13_WoL_index_filtered',\n",
+ " 'post_left_sub8_S14_WoL_index_filtered',\n",
+ " 'post_left_sub9_S110_WoL_index_filtered'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['post_left_sub11_S16_WoL_index_filtered',\n",
+ " 'post_left_sub5_S12_WoL_index_filtered',\n",
+ " 'post_left_sub6_S13_WoL_index_filtered',\n",
+ " 'post_left_sub8_S14_WoL_index_filtered',\n",
+ " 'post_left_sub9_S110_WoL_index_filtered'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['post_left_sub11_S16_WoL_index_filtered',\n",
+ " 'post_left_sub5_S12_WoL_index_filtered',\n",
+ " 'post_left_sub6_S13_WoL_index_filtered',\n",
+ " 'post_left_sub8_S14_WoL_index_filtered',\n",
+ " 'post_left_sub9_S110_WoL_index_filtered'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['post_left_sub11_S16_WoL_index_filtered',\n",
+ " 'post_left_sub5_S12_WoL_index_filtered',\n",
+ " 'post_left_sub6_S13_WoL_index_filtered',\n",
+ " 'post_left_sub8_S14_WoL_index_filtered',\n",
+ " 'post_left_sub9_S110_WoL_index_filtered'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['post_left_sub11_S16_WoL_index_filtered',\n",
+ " 'post_left_sub5_S12_WoL_index_filtered',\n",
+ " 'post_left_sub6_S13_WoL_index_filtered',\n",
+ " 'post_left_sub8_S14_WoL_index_filtered',\n",
+ " 'post_left_sub9_S110_WoL_index_filtered'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['post_left_sub11_S16_WoL_index_filtered',\n",
+ " 'post_left_sub5_S12_WoL_index_filtered',\n",
+ " 'post_left_sub6_S13_WoL_index_filtered',\n",
+ " 'post_left_sub8_S14_WoL_index_filtered',\n",
+ " 'post_left_sub9_S110_WoL_index_filtered'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['post_left_sub11_S16_WoL_index_filtered',\n",
+ " 'post_left_sub5_S12_WoL_index_filtered',\n",
+ " 'post_left_sub6_S13_WoL_index_filtered',\n",
+ " 'post_left_sub8_S14_WoL_index_filtered',\n",
+ " 'post_left_sub9_S110_WoL_index_filtered'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['post_left_sub11_S16_WoL_index_filtered',\n",
+ " 'post_left_sub5_S12_WoL_index_filtered',\n",
+ " 'post_left_sub6_S13_WoL_index_filtered',\n",
+ " 'post_left_sub8_S14_WoL_index_filtered',\n",
+ " 'post_left_sub9_S110_WoL_index_filtered'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['post_left_sub11_S16_WoL_index_filtered',\n",
+ " 'post_left_sub5_S12_WoL_index_filtered',\n",
+ " 'post_left_sub6_S13_WoL_index_filtered',\n",
+ " 'post_left_sub8_S14_WoL_index_filtered',\n",
+ " 'post_left_sub9_S110_WoL_index_filtered'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['post_left_sub11_S16_WoL_index_filtered',\n",
+ " 'post_left_sub5_S12_WoL_index_filtered',\n",
+ " 'post_left_sub6_S13_WoL_index_filtered',\n",
+ " 'post_left_sub8_S14_WoL_index_filtered',\n",
+ " 'post_left_sub9_S110_WoL_index_filtered'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['post_left_sub11_S16_WoL_index_filtered',\n",
+ " 'post_left_sub5_S12_WoL_index_filtered',\n",
+ " 'post_left_sub6_S13_WoL_index_filtered',\n",
+ " 'post_left_sub8_S14_WoL_index_filtered',\n",
+ " 'post_left_sub9_S110_WoL_index_filtered'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['post_left_sub11_S16_WoL_index_filtered',\n",
+ " 'post_left_sub5_S12_WoL_index_filtered',\n",
+ " 'post_left_sub6_S13_WoL_index_filtered',\n",
+ " 'post_left_sub8_S14_WoL_index_filtered',\n",
+ " 'post_left_sub9_S110_WoL_index_filtered'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['post_left_sub11_S16_WoL_index_filtered',\n",
+ " 'post_left_sub5_S12_WoL_index_filtered',\n",
+ " 'post_left_sub6_S13_WoL_index_filtered',\n",
+ " 'post_left_sub8_S14_WoL_index_filtered',\n",
+ " 'post_left_sub9_S110_WoL_index_filtered'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['post_left_sub11_S16_WoL_index_filtered',\n",
+ " 'post_left_sub5_S12_WoL_index_filtered',\n",
+ " 'post_left_sub6_S13_WoL_index_filtered',\n",
+ " 'post_left_sub8_S14_WoL_index_filtered',\n",
+ " 'post_left_sub9_S110_WoL_index_filtered'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['post_left_sub11_S16_WoL_index_filtered',\n",
+ " 'post_left_sub5_S12_WoL_index_filtered',\n",
+ " 'post_left_sub6_S13_WoL_index_filtered',\n",
+ " 'post_left_sub8_S14_WoL_index_filtered',\n",
+ " 'post_left_sub9_S110_WoL_index_filtered'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['post_left_sub11_S16_WoL_index_filtered',\n",
+ " 'post_left_sub5_S12_WoL_index_filtered',\n",
+ " 'post_left_sub6_S13_WoL_index_filtered',\n",
+ " 'post_left_sub8_S14_WoL_index_filtered',\n",
+ " 'post_left_sub9_S110_WoL_index_filtered'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['post_left_sub11_S16_WoL_index_filtered',\n",
+ " 'post_left_sub5_S12_WoL_index_filtered',\n",
+ " 'post_left_sub6_S13_WoL_index_filtered',\n",
+ " 'post_left_sub8_S14_WoL_index_filtered',\n",
+ " 'post_left_sub9_S110_WoL_index_filtered'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['post_left_sub11_S16_WoL_index_filtered',\n",
+ " 'post_left_sub5_S12_WoL_index_filtered',\n",
+ " 'post_left_sub6_S13_WoL_index_filtered',\n",
+ " 'post_left_sub8_S14_WoL_index_filtered',\n",
+ " 'post_left_sub9_S110_WoL_index_filtered'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['post_left_sub11_S16_WoL_index_filtered',\n",
+ " 'post_left_sub5_S12_WoL_index_filtered',\n",
+ " 'post_left_sub6_S13_WoL_index_filtered',\n",
+ " 'post_left_sub8_S14_WoL_index_filtered',\n",
+ " 'post_left_sub9_S110_WoL_index_filtered'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['post_left_sub11_S16_WoL_index_filtered',\n",
+ " 'post_left_sub5_S12_WoL_index_filtered',\n",
+ " 'post_left_sub6_S13_WoL_index_filtered',\n",
+ " 'post_left_sub8_S14_WoL_index_filtered',\n",
+ " 'post_left_sub9_S110_WoL_index_filtered'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['post_left_sub11_S16_WoL_index_filtered',\n",
+ " 'post_left_sub5_S12_WoL_index_filtered',\n",
+ " 'post_left_sub6_S13_WoL_index_filtered',\n",
+ " 'post_left_sub8_S14_WoL_index_filtered',\n",
+ " 'post_left_sub9_S110_WoL_index_filtered'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['post_left_sub11_S16_WoL_index_filtered',\n",
+ " 'post_left_sub5_S12_WoL_index_filtered',\n",
+ " 'post_left_sub6_S13_WoL_index_filtered',\n",
+ " 'post_left_sub8_S14_WoL_index_filtered',\n",
+ " 'post_left_sub9_S110_WoL_index_filtered'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['post_left_sub11_S16_WoL_index_filtered',\n",
+ " 'post_left_sub5_S12_WoL_index_filtered',\n",
+ " 'post_left_sub6_S13_WoL_index_filtered',\n",
+ " 'post_left_sub8_S14_WoL_index_filtered',\n",
+ " 'post_left_sub9_S110_WoL_index_filtered'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['post_left_sub11_S16_WoL_index_filtered',\n",
+ " 'post_left_sub5_S12_WoL_index_filtered',\n",
+ " 'post_left_sub6_S13_WoL_index_filtered',\n",
+ " 'post_left_sub8_S14_WoL_index_filtered',\n",
+ " 'post_left_sub9_S110_WoL_index_filtered'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['post_left_sub11_S16_WoL_index_filtered',\n",
+ " 'post_left_sub5_S12_WoL_index_filtered',\n",
+ " 'post_left_sub6_S13_WoL_index_filtered',\n",
+ " 'post_left_sub8_S14_WoL_index_filtered',\n",
+ " 'post_left_sub9_S110_WoL_index_filtered'],\n",
+ " dtype='object', name='sampleid')\n"
+ ]
+ },
+ {
+ "data": {
+ "text/html": [
+ "\n",
+ "\n",
+ "
\n",
+ " \n",
+ " \n",
+ " | \n",
+ " 24hControl_vs_24hChemical_tstat | \n",
+ " 24hControl__vs_24hChemical_pvalue | \n",
+ " 24hControl_vs_24hChemical_upperCI | \n",
+ " 24hControl_vs_24hChemical_lowerCI | \n",
+ "
\n",
+ " \n",
+ " \n",
+ " \n",
+ " | log(C.acnes/C.afermentans) | \n",
+ " -0.346499 | \n",
+ " 0.746428 | \n",
+ " -0.194093 | \n",
+ " 0.151023 | \n",
+ "
\n",
+ " \n",
+ " | log(M.luteus/C.afermentans) | \n",
+ " -0.076048 | \n",
+ " 0.943033 | \n",
+ " -0.463871 | \n",
+ " 0.439137 | \n",
+ "
\n",
+ " \n",
+ " | log(S.aureus/C.afermentans) | \n",
+ " 0.497989 | \n",
+ " 0.644627 | \n",
+ " -0.174774 | \n",
+ " 0.251172 | \n",
+ "
\n",
+ " \n",
+ " | log(S.capitis/C.afermentans) | \n",
+ " 0.822257 | \n",
+ " 0.457119 | \n",
+ " -0.438598 | \n",
+ " 0.807692 | \n",
+ "
\n",
+ " \n",
+ " | log(S.epidermidis/C.afermentans) | \n",
+ " 0.928363 | \n",
+ " 0.405768 | \n",
+ " -0.307996 | \n",
+ " 0.617433 | \n",
+ "
\n",
+ " \n",
+ " | log(S.hominis/C.afermentans) | \n",
+ " -0.700531 | \n",
+ " 0.522202 | \n",
+ " -1.043407 | \n",
+ " 0.622962 | \n",
+ "
\n",
+ " \n",
+ " | log(S.warneri/C.afermentans) | \n",
+ " 1.669434 | \n",
+ " 0.170354 | \n",
+ " -0.110073 | \n",
+ " 0.442065 | \n",
+ "
\n",
+ " \n",
+ " | log(S.mitis/C.afermentans) | \n",
+ " -1.047233 | \n",
+ " 0.354098 | \n",
+ " -1.786511 | \n",
+ " 0.807928 | \n",
+ "
\n",
+ " \n",
+ "
\n",
+ "
"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "# table normalization by column\n",
+ "table = table.loc[:, table.sum(axis=0)>0]\n",
+ "table = table.apply(lambda x: x / x.sum(), axis=1)\n",
+ "\n",
+ "reference_species = 'Corynebacterium afermentans'\n",
+ "\n",
+ "def balance_f(x):\n",
+ " lr = [np.log(x['Cutibacterium acnes']+0.001) - np.log(x[reference_species]+0.001),\n",
+ " np.log(x['Micrococcus luteus']+0.001) - np.log(x[reference_species]+0.001),\n",
+ " np.log(x['Staphylococcus aureus']+0.001) - np.log(x[reference_species]+0.001),\n",
+ " np.log(x['Staphylococcus capitis']+0.001) - np.log(x[reference_species]+0.001),\n",
+ " np.log(x['Staphylococcus epidermidis']+0.001) - np.log(x[reference_species]+0.001),\n",
+ " np.log(x['Staphylococcus hominis']+0.001) - np.log(x[reference_species]+0.001),\n",
+ " np.log(x['Staphylococcus warneri']+0.001) - np.log(x[reference_species]+0.001),\n",
+ " np.log(x['Streptococcus mitis']+0.001) - np.log(x[reference_species]+0.001)\n",
+ " ]\n",
+ " cols = [\n",
+ " 'log(C.acnes/C.afermentans)',\n",
+ " 'log(M.luteus/C.afermentans)',\n",
+ " 'log(S.aureus/C.afermentans)',\n",
+ " 'log(S.capitis/C.afermentans)',\n",
+ " 'log(S.epidermidis/C.afermentans)',\n",
+ " 'log(S.hominis/C.afermentans)',\n",
+ " 'log(S.warneri/C.afermentans)',\n",
+ " 'log(S.mitis/C.afermentans)'\n",
+ " ]\n",
+ " return pd.Series(lr, index=cols)\n",
+ "\n",
+ "balances = table.apply(balance_f, axis = 1)\n",
+ "\n",
+ "f = lambda x: paired_t_test(category='Control', md=metadata, df=x)[0]\n",
+ "tstats = balances.apply(f, axis=0)\n",
+ "\n",
+ "f = lambda x: paired_t_test(category='Control', md=metadata, df=x)[1]\n",
+ "tpvals = balances.apply(f, axis=0)\n",
+ "\n",
+ "f = lambda x: paired_t_test(category='Control', md=metadata, df=x)[2]\n",
+ "tlower = balances.apply(f, axis=0)\n",
+ "\n",
+ "f = lambda x: paired_t_test(category='Control', md=metadata, df=x)[3]\n",
+ "tupper = balances.apply(f, axis=0)\n",
+ "\n",
+ "stats = pd.DataFrame(\n",
+ " {\n",
+ " '24hControl_vs_24hChemical_tstat': tstats,\n",
+ " '24hControl__vs_24hChemical_pvalue': tpvals,\n",
+ " '24hControl_vs_24hChemical_upperCI': tlower,\n",
+ " '24hControl_vs_24hChemical_lowerCI': tupper\n",
+ " }\n",
+ ")\n",
+ "\n",
+ "data = pd.merge(balances, metadata, left_index=True, right_index=True)\n",
+ "\n",
+ "subdata = data[[\n",
+ " 'log(C.acnes/C.afermentans)',\n",
+ " 'log(M.luteus/C.afermentans)',\n",
+ " 'log(S.aureus/C.afermentans)',\n",
+ " 'log(S.capitis/C.afermentans)',\n",
+ " 'log(S.epidermidis/C.afermentans)',\n",
+ " 'log(S.hominis/C.afermentans)',\n",
+ " 'log(S.warneri/C.afermentans)',\n",
+ " 'log(S.mitis/C.afermentans)',\n",
+ " 'Control']]\n",
+ "\n",
+ "df = pd.melt(subdata, id_vars = 'Control', var_name = 'balance')\n",
+ "\n",
+ "# get the palette color name\n",
+ "deep_palette = sns.color_palette(\"Set1\")\n",
+ "\n",
+ "# Extract the first two colors\n",
+ "color1 = deep_palette[0]\n",
+ "color2 = deep_palette[1]\n",
+ "\n",
+ "#, notch=True\n",
+ "sns.set(style=\"white\", font_scale=1.5, font=\"sans-serif\", palette=\"Set1\")\n",
+ "\n",
+ "# get a label for outlier,figure out why there is error bar \n",
+ "# customize which counts as outlier\n",
+ "g = sns.catplot(x='value', y='balance', data=df, hue='Control', kind='box',\n",
+ " hue_order=['left', 'right'], n_boot = 1000,\n",
+ " height = 10, aspect = 1.8, legend = False, whis = [0,100])\n",
+ "g.set_xlabels('log fold change', fontsize=14)\n",
+ "g.set_ylabels('', fontsize=14)\n",
+ "plt.legend(fontsize = 20, frameon=False)\n",
+ "# plt.legend(loc='lower right', labels=['Control', 'Treatment'])\n",
+ "plt.legend(fontsize=20, frameon=False, loc='lower right')\n",
+ "plt.xticks(fontsize=12)\n",
+ "\n",
+ "ax = plt.gca()\n",
+ "leg = ax.legend(fontsize=20, frameon=False, loc='lower right')\n",
+ "leg.legendHandles[0].set_linewidth(8)\n",
+ "leg.legendHandles[0].set_color(color1)\n",
+ "leg.legendHandles[1].set_linewidth(8)\n",
+ "leg.legendHandles[1].set_color(color2)\n",
+ "leg.get_frame().set_linewidth(0)\n",
+ "#plt.xticks(rotation=60)\n",
+ "\n",
+ "# g.savefig('./out_sles/fig_sles.svg', format='svg')\n",
+ "# stats.to_csv('./out_c_afermentans/stats_day7_rl.csv')\n",
+ "stats"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "95f2acab-fbed-4ba2-9ff2-87334e09af9d",
+ "metadata": {},
+ "outputs": [],
+ "source": []
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "qiime2-2020.6-system",
+ "language": "python",
+ "name": "qiime2-2020.6-system"
+ },
+ "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.10"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/figure_generation_invitro_afermentans_reordered.ipynb b/figure_generation_invitro_afermentans_reordered.ipynb
new file mode 100644
index 0000000..192774d
--- /dev/null
+++ b/figure_generation_invitro_afermentans_reordered.ipynb
@@ -0,0 +1,1124 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": 39,
+ "id": "7998a5a1-e3de-42b7-980c-626c1ae24168",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# ----------------------------------------------------------------------------\n",
+ "# SkinCom Analysis - Perform t-test & boxplots for log-ratios\n",
+ "# Author: Sherlyn Weng\n",
+ "# Version: 1.0.0\n",
+ "# Mmaintainer: Sherlyn Weng\n",
+ "# Email: y1weng@ucsd.edu\n",
+ "# This code is adapted from https://github.com/knightlab-analyses/reference-frames\n",
+ "# ----------------------------------------------------------------------------"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 40,
+ "id": "005db807-d8b9-4eb4-893b-ec07aba48ec3",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import numpy as np\n",
+ "import pandas as pd\n",
+ "import seaborn as sns\n",
+ "import matplotlib.pyplot as plt\n",
+ "import matplotlib as mpl\n",
+ "from biom import load_table\n",
+ "from scipy.stats import ttest_rel, wilcoxon, ttest_ind, pearsonr, spearmanr\n",
+ "from skbio.stats.composition import multiplicative_replacement\n",
+ "import matplotlib.gridspec as gridspec\n",
+ "\n",
+ "from collections import namedtuple\n",
+ "\n",
+ "%matplotlib inline\n",
+ "\n",
+ "import warnings\n",
+ "warnings.filterwarnings('ignore')\n",
+ "\n",
+ "\"\"\"\n",
+ "Here, we have implemented an extended version of scipy's paired t-test function.\n",
+ "see https://github.com/scipy/scipy/blob/master/scipy/stats/mstats_basic.py#L1085\n",
+ "\n",
+ "This should really be pushed to scipy -- any volunteers?\n",
+ "\"\"\"\n",
+ "import numpy as np\n",
+ "from scipy.stats import t\n",
+ "from scipy.special import betainc\n",
+ "from collections import namedtuple\n",
+ "\n",
+ "\n",
+ "Ttest_relResult = namedtuple('Ttest_relResult', \n",
+ " ('statistic', 'pvalue', 'lower_CI', 'upper_CI', 'ddof'))\n",
+ "\n",
+ "\n",
+ "def ttest_rel(x, y):\n",
+ " \n",
+ " D = x - y\n",
+ " n = len(D)\n",
+ " \n",
+ " sd = np.std(D, ddof=1) \n",
+ " se = sd / np.sqrt(n)\n",
+ " mu = np.mean(D)\n",
+ " T = mu / se\n",
+ " ddof = n-1\n",
+ " ci = se * t.ppf( 0.05 / 2, ddof) \n",
+ " lower = mu - np.abs(ci)\n",
+ " upper = mu + np.abs(ci)\n",
+ " \n",
+ " prob = betainc(0.5*ddof, 0.5, ddof/(ddof + T*T))\n",
+ " return Ttest_relResult(T, prob, lower, upper, int(ddof))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 41,
+ "id": "2e265a3f-520c-43bc-955b-20b179fe1cee",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def paired_t_test(category, md, df):\n",
+ " \"\"\" Performs paired t-test on a dataframe with complementing\n",
+ " sample metadata.\n",
+ "\n",
+ " Parameters\n",
+ " ----------\n",
+ " category : str\n",
+ " Category to test in the dataframe df\n",
+ " md : pd.DataFrame\n",
+ " Sample metadata specifying groupings\n",
+ " df : pd.Dataframe\n",
+ " Counts of microbes across all samples\n",
+ " \n",
+ " Returns\n",
+ " -------\n",
+ " T : float\n",
+ " T-statistic\n",
+ " prob : float\n",
+ " p-value\n",
+ " lower : float\n",
+ " lower 95% CI bound\n",
+ " upper : float\n",
+ " upper 95% CI bound\n",
+ " ddof : int\n",
+ " degrees of freedom\n",
+ " \"\"\"\n",
+ " cats = list(md[category].value_counts().index)\n",
+ " src = md[md[category] == cats[0]].index\n",
+ " print('hellooo')\n",
+ " print(src)\n",
+ " dest = md[md[category] == cats[1]].index\n",
+ " return ttest_rel(df.loc[dest].values, df.loc[src].values)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "5242c6a2-b6b2-497e-aa34-225f3a9c1cab",
+ "metadata": {},
+ "source": [
+ "## Input Directory"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 42,
+ "id": "0f656bb7-fbf6-4b0e-909b-d94835646df4",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "input_dir = \"../20231128_run_songbird_skincom\""
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "f9b78015-1a9a-4d8e-97d0-e7ff85a1e82c",
+ "metadata": {},
+ "source": [
+ "## Variables"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "1164688b-1357-4d73-bedb-66dc5b71386a",
+ "metadata": {
+ "tags": []
+ },
+ "source": [
+ "## SLES"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 43,
+ "id": "d993c044-e503-4682-9279-dbde03602d26",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "hellooo\n",
+ "Index(['Comm-1_S4_L001', 'Comm-1_S4_L002', 'Comm-2_S14_L001',\n",
+ " 'Comm-2_S14_L002', 'Comm-3_S23_L001', 'Comm-3_S23_L002',\n",
+ " 'Comm-4_S32_L001', 'Comm-4_S32_L002'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hellooo\n",
+ "Index(['Comm-1_S4_L001', 'Comm-1_S4_L002', 'Comm-2_S14_L001',\n",
+ " 'Comm-2_S14_L002', 'Comm-3_S23_L001', 'Comm-3_S23_L002',\n",
+ " 'Comm-4_S32_L001', 'Comm-4_S32_L002'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hellooo\n",
+ "Index(['Comm-1_S4_L001', 'Comm-1_S4_L002', 'Comm-2_S14_L001',\n",
+ " 'Comm-2_S14_L002', 'Comm-3_S23_L001', 'Comm-3_S23_L002',\n",
+ " 'Comm-4_S32_L001', 'Comm-4_S32_L002'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hellooo\n",
+ "Index(['Comm-1_S4_L001', 'Comm-1_S4_L002', 'Comm-2_S14_L001',\n",
+ " 'Comm-2_S14_L002', 'Comm-3_S23_L001', 'Comm-3_S23_L002',\n",
+ " 'Comm-4_S32_L001', 'Comm-4_S32_L002'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hellooo\n",
+ "Index(['Comm-1_S4_L001', 'Comm-1_S4_L002', 'Comm-2_S14_L001',\n",
+ " 'Comm-2_S14_L002', 'Comm-3_S23_L001', 'Comm-3_S23_L002',\n",
+ " 'Comm-4_S32_L001', 'Comm-4_S32_L002'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hellooo\n",
+ "Index(['Comm-1_S4_L001', 'Comm-1_S4_L002', 'Comm-2_S14_L001',\n",
+ " 'Comm-2_S14_L002', 'Comm-3_S23_L001', 'Comm-3_S23_L002',\n",
+ " 'Comm-4_S32_L001', 'Comm-4_S32_L002'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hellooo\n",
+ "Index(['Comm-1_S4_L001', 'Comm-1_S4_L002', 'Comm-2_S14_L001',\n",
+ " 'Comm-2_S14_L002', 'Comm-3_S23_L001', 'Comm-3_S23_L002',\n",
+ " 'Comm-4_S32_L001', 'Comm-4_S32_L002'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hellooo\n",
+ "Index(['Comm-1_S4_L001', 'Comm-1_S4_L002', 'Comm-2_S14_L001',\n",
+ " 'Comm-2_S14_L002', 'Comm-3_S23_L001', 'Comm-3_S23_L002',\n",
+ " 'Comm-4_S32_L001', 'Comm-4_S32_L002'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hellooo\n",
+ "Index(['Comm-1_S4_L001', 'Comm-1_S4_L002', 'Comm-2_S14_L001',\n",
+ " 'Comm-2_S14_L002', 'Comm-3_S23_L001', 'Comm-3_S23_L002',\n",
+ " 'Comm-4_S32_L001', 'Comm-4_S32_L002'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hellooo\n",
+ "Index(['Comm-1_S4_L001', 'Comm-1_S4_L002', 'Comm-2_S14_L001',\n",
+ " 'Comm-2_S14_L002', 'Comm-3_S23_L001', 'Comm-3_S23_L002',\n",
+ " 'Comm-4_S32_L001', 'Comm-4_S32_L002'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hellooo\n",
+ "Index(['Comm-1_S4_L001', 'Comm-1_S4_L002', 'Comm-2_S14_L001',\n",
+ " 'Comm-2_S14_L002', 'Comm-3_S23_L001', 'Comm-3_S23_L002',\n",
+ " 'Comm-4_S32_L001', 'Comm-4_S32_L002'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hellooo\n",
+ "Index(['Comm-1_S4_L001', 'Comm-1_S4_L002', 'Comm-2_S14_L001',\n",
+ " 'Comm-2_S14_L002', 'Comm-3_S23_L001', 'Comm-3_S23_L002',\n",
+ " 'Comm-4_S32_L001', 'Comm-4_S32_L002'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hellooo\n",
+ "Index(['Comm-1_S4_L001', 'Comm-1_S4_L002', 'Comm-2_S14_L001',\n",
+ " 'Comm-2_S14_L002', 'Comm-3_S23_L001', 'Comm-3_S23_L002',\n",
+ " 'Comm-4_S32_L001', 'Comm-4_S32_L002'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hellooo\n",
+ "Index(['Comm-1_S4_L001', 'Comm-1_S4_L002', 'Comm-2_S14_L001',\n",
+ " 'Comm-2_S14_L002', 'Comm-3_S23_L001', 'Comm-3_S23_L002',\n",
+ " 'Comm-4_S32_L001', 'Comm-4_S32_L002'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hellooo\n",
+ "Index(['Comm-1_S4_L001', 'Comm-1_S4_L002', 'Comm-2_S14_L001',\n",
+ " 'Comm-2_S14_L002', 'Comm-3_S23_L001', 'Comm-3_S23_L002',\n",
+ " 'Comm-4_S32_L001', 'Comm-4_S32_L002'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hellooo\n",
+ "Index(['Comm-1_S4_L001', 'Comm-1_S4_L002', 'Comm-2_S14_L001',\n",
+ " 'Comm-2_S14_L002', 'Comm-3_S23_L001', 'Comm-3_S23_L002',\n",
+ " 'Comm-4_S32_L001', 'Comm-4_S32_L002'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hellooo\n",
+ "Index(['Comm-1_S4_L001', 'Comm-1_S4_L002', 'Comm-2_S14_L001',\n",
+ " 'Comm-2_S14_L002', 'Comm-3_S23_L001', 'Comm-3_S23_L002',\n",
+ " 'Comm-4_S32_L001', 'Comm-4_S32_L002'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hellooo\n",
+ "Index(['Comm-1_S4_L001', 'Comm-1_S4_L002', 'Comm-2_S14_L001',\n",
+ " 'Comm-2_S14_L002', 'Comm-3_S23_L001', 'Comm-3_S23_L002',\n",
+ " 'Comm-4_S32_L001', 'Comm-4_S32_L002'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hellooo\n",
+ "Index(['Comm-1_S4_L001', 'Comm-1_S4_L002', 'Comm-2_S14_L001',\n",
+ " 'Comm-2_S14_L002', 'Comm-3_S23_L001', 'Comm-3_S23_L002',\n",
+ " 'Comm-4_S32_L001', 'Comm-4_S32_L002'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hellooo\n",
+ "Index(['Comm-1_S4_L001', 'Comm-1_S4_L002', 'Comm-2_S14_L001',\n",
+ " 'Comm-2_S14_L002', 'Comm-3_S23_L001', 'Comm-3_S23_L002',\n",
+ " 'Comm-4_S32_L001', 'Comm-4_S32_L002'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hellooo\n",
+ "Index(['Comm-1_S4_L001', 'Comm-1_S4_L002', 'Comm-2_S14_L001',\n",
+ " 'Comm-2_S14_L002', 'Comm-3_S23_L001', 'Comm-3_S23_L002',\n",
+ " 'Comm-4_S32_L001', 'Comm-4_S32_L002'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hellooo\n",
+ "Index(['Comm-1_S4_L001', 'Comm-1_S4_L002', 'Comm-2_S14_L001',\n",
+ " 'Comm-2_S14_L002', 'Comm-3_S23_L001', 'Comm-3_S23_L002',\n",
+ " 'Comm-4_S32_L001', 'Comm-4_S32_L002'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hellooo\n",
+ "Index(['Comm-1_S4_L001', 'Comm-1_S4_L002', 'Comm-2_S14_L001',\n",
+ " 'Comm-2_S14_L002', 'Comm-3_S23_L001', 'Comm-3_S23_L002',\n",
+ " 'Comm-4_S32_L001', 'Comm-4_S32_L002'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hellooo\n",
+ "Index(['Comm-1_S4_L001', 'Comm-1_S4_L002', 'Comm-2_S14_L001',\n",
+ " 'Comm-2_S14_L002', 'Comm-3_S23_L001', 'Comm-3_S23_L002',\n",
+ " 'Comm-4_S32_L001', 'Comm-4_S32_L002'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hellooo\n",
+ "Index(['Comm-1_S4_L001', 'Comm-1_S4_L002', 'Comm-2_S14_L001',\n",
+ " 'Comm-2_S14_L002', 'Comm-3_S23_L001', 'Comm-3_S23_L002',\n",
+ " 'Comm-4_S32_L001', 'Comm-4_S32_L002'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hellooo\n",
+ "Index(['Comm-1_S4_L001', 'Comm-1_S4_L002', 'Comm-2_S14_L001',\n",
+ " 'Comm-2_S14_L002', 'Comm-3_S23_L001', 'Comm-3_S23_L002',\n",
+ " 'Comm-4_S32_L001', 'Comm-4_S32_L002'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hellooo\n",
+ "Index(['Comm-1_S4_L001', 'Comm-1_S4_L002', 'Comm-2_S14_L001',\n",
+ " 'Comm-2_S14_L002', 'Comm-3_S23_L001', 'Comm-3_S23_L002',\n",
+ " 'Comm-4_S32_L001', 'Comm-4_S32_L002'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hellooo\n",
+ "Index(['Comm-1_S4_L001', 'Comm-1_S4_L002', 'Comm-2_S14_L001',\n",
+ " 'Comm-2_S14_L002', 'Comm-3_S23_L001', 'Comm-3_S23_L002',\n",
+ " 'Comm-4_S32_L001', 'Comm-4_S32_L002'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hellooo\n",
+ "Index(['Comm-1_S4_L001', 'Comm-1_S4_L002', 'Comm-2_S14_L001',\n",
+ " 'Comm-2_S14_L002', 'Comm-3_S23_L001', 'Comm-3_S23_L002',\n",
+ " 'Comm-4_S32_L001', 'Comm-4_S32_L002'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hellooo\n",
+ "Index(['Comm-1_S4_L001', 'Comm-1_S4_L002', 'Comm-2_S14_L001',\n",
+ " 'Comm-2_S14_L002', 'Comm-3_S23_L001', 'Comm-3_S23_L002',\n",
+ " 'Comm-4_S32_L001', 'Comm-4_S32_L002'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hellooo\n",
+ "Index(['Comm-1_S4_L001', 'Comm-1_S4_L002', 'Comm-2_S14_L001',\n",
+ " 'Comm-2_S14_L002', 'Comm-3_S23_L001', 'Comm-3_S23_L002',\n",
+ " 'Comm-4_S32_L001', 'Comm-4_S32_L002'],\n",
+ " dtype='object', name='sampleid')\n",
+ "hellooo\n",
+ "Index(['Comm-1_S4_L001', 'Comm-1_S4_L002', 'Comm-2_S14_L001',\n",
+ " 'Comm-2_S14_L002', 'Comm-3_S23_L001', 'Comm-3_S23_L002',\n",
+ " 'Comm-4_S32_L001', 'Comm-4_S32_L002'],\n",
+ " dtype='object', name='sampleid')\n"
+ ]
+ },
+ {
+ "data": {
+ "text/html": [
+ "\n",
+ "\n",
+ "
\n",
+ " \n",
+ " \n",
+ " | \n",
+ " 24hControl_vs_24hChemical_tstat | \n",
+ " 24hControl__vs_24hChemical_pvalue | \n",
+ " 24hControl_vs_24hChemical_upperCI | \n",
+ " 24hControl_vs_24hChemical_lowerCI | \n",
+ "
\n",
+ " \n",
+ " \n",
+ " \n",
+ " | log(C.acnes/C.afermentans) | \n",
+ " -4.546796 | \n",
+ " 0.002646 | \n",
+ " -3.477390 | \n",
+ " -1.097931 | \n",
+ "
\n",
+ " \n",
+ " | log(M.luteus/C.afermentans) | \n",
+ " 1.786337 | \n",
+ " 0.117204 | \n",
+ " -0.013707 | \n",
+ " 0.098391 | \n",
+ "
\n",
+ " \n",
+ " | log(S.aureus/C.afermentans) | \n",
+ " -12.891296 | \n",
+ " 0.000004 | \n",
+ " -4.310889 | \n",
+ " -2.974538 | \n",
+ "
\n",
+ " \n",
+ " | log(S.capitis/C.afermentans) | \n",
+ " 0.706425 | \n",
+ " 0.502753 | \n",
+ " -0.041647 | \n",
+ " 0.077132 | \n",
+ "
\n",
+ " \n",
+ " | log(S.epidermidis/C.afermentans) | \n",
+ " 1.257191 | \n",
+ " 0.249005 | \n",
+ " -0.074636 | \n",
+ " 0.244093 | \n",
+ "
\n",
+ " \n",
+ " | log(S.hominis/C.afermentans) | \n",
+ " 9.084779 | \n",
+ " 0.000040 | \n",
+ " 0.452786 | \n",
+ " 0.771429 | \n",
+ "
\n",
+ " \n",
+ " | log(S.warneri/C.afermentans) | \n",
+ " 8.768258 | \n",
+ " 0.000051 | \n",
+ " 0.333996 | \n",
+ " 0.580661 | \n",
+ "
\n",
+ " \n",
+ " | log(S.mitis/C.afermentans) | \n",
+ " -7.348948 | \n",
+ " 0.000156 | \n",
+ " -0.141701 | \n",
+ " -0.072711 | \n",
+ "
\n",
+ " \n",
+ "
\n",
+ "
"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "# load md\n",
+ "metadata = pd.read_table(f'{input_dir}/data/syncom_sles_md.tsv', index_col = 0)\n",
+ "\n",
+ "#load samples\n",
+ "table = load_table(f'{input_dir}/data/syncom_sles.biom').to_dataframe().T\n",
+ "\n",
+ "# load multinomial ranks (differentials)\n",
+ "beta = pd.read_csv(f'{input_dir}/out_sles/differentials.tsv', sep='\\t', index_col = 0)\n",
+ "beta = beta.iloc[1:]\n",
+ "\n",
+ "# table normalization by column\n",
+ "table = table.loc[:, table.sum(axis=0)>0]\n",
+ "table = table.apply(lambda x: x / x.sum(), axis=1)\n",
+ "\n",
+ "reference_species = 'Corynebacterium afermentans'\n",
+ "\n",
+ "def balance_f(x):\n",
+ " lr = [np.log(x['Cutibacterium acnes']+0.001) - np.log(x[reference_species]+0.001),\n",
+ " np.log(x['Micrococcus luteus']+0.001) - np.log(x[reference_species]+0.001),\n",
+ " np.log(x['Staphylococcus aureus']+0.001) - np.log(x[reference_species]+0.001),\n",
+ " np.log(x['Staphylococcus capitis']+0.001) - np.log(x[reference_species]+0.001),\n",
+ " np.log(x['Staphylococcus epidermidis']+0.001) - np.log(x[reference_species]+0.001),\n",
+ " np.log(x['Staphylococcus hominis']+0.001) - np.log(x[reference_species]+0.001),\n",
+ " np.log(x['Staphylococcus warneri']+0.001) - np.log(x[reference_species]+0.001),\n",
+ " np.log(x['Streptococcus mitis']+0.001) - np.log(x[reference_species]+0.001)\n",
+ " ]\n",
+ " cols = [\n",
+ " 'log(C.acnes/C.afermentans)',\n",
+ " 'log(M.luteus/C.afermentans)',\n",
+ " 'log(S.aureus/C.afermentans)',\n",
+ " 'log(S.capitis/C.afermentans)',\n",
+ " 'log(S.epidermidis/C.afermentans)',\n",
+ " 'log(S.hominis/C.afermentans)',\n",
+ " 'log(S.warneri/C.afermentans)',\n",
+ " 'log(S.mitis/C.afermentans)'\n",
+ " ]\n",
+ " return pd.Series(lr, index=cols)\n",
+ "\n",
+ "balances = table.apply(balance_f, axis = 1)\n",
+ "\n",
+ "f = lambda x: paired_t_test(category='treatment', md=metadata, df=x)[0]\n",
+ "tstats = balances.apply(f, axis=0)\n",
+ "\n",
+ "f = lambda x: paired_t_test(category='treatment', md=metadata, df=x)[1]\n",
+ "tpvals = balances.apply(f, axis=0)\n",
+ "\n",
+ "f = lambda x: paired_t_test(category='treatment', md=metadata, df=x)[2]\n",
+ "tlower = balances.apply(f, axis=0)\n",
+ "\n",
+ "f = lambda x: paired_t_test(category='treatment', md=metadata, df=x)[3]\n",
+ "tupper = balances.apply(f, axis=0)\n",
+ "\n",
+ "stats = pd.DataFrame(\n",
+ " {\n",
+ " '24hControl_vs_24hChemical_tstat': tstats,\n",
+ " '24hControl__vs_24hChemical_pvalue': tpvals,\n",
+ " '24hControl_vs_24hChemical_upperCI': tlower,\n",
+ " '24hControl_vs_24hChemical_lowerCI': tupper\n",
+ " }\n",
+ ")\n",
+ "\n",
+ "data = pd.merge(balances, metadata, left_index=True, right_index=True)\n",
+ "\n",
+ "subdata = data[[\n",
+ " 'log(C.acnes/C.afermentans)',\n",
+ " 'log(M.luteus/C.afermentans)',\n",
+ " 'log(S.aureus/C.afermentans)',\n",
+ " 'log(S.capitis/C.afermentans)',\n",
+ " 'log(S.epidermidis/C.afermentans)',\n",
+ " 'log(S.hominis/C.afermentans)',\n",
+ " 'log(S.warneri/C.afermentans)',\n",
+ " 'log(S.mitis/C.afermentans)',\n",
+ " 'treatment']]\n",
+ "\n",
+ "df = pd.melt(subdata, id_vars = 'treatment', var_name = 'balance')\n",
+ "\n",
+ "# get the palette color name\n",
+ "deep_palette = sns.color_palette(\"Set1\")\n",
+ "\n",
+ "# Extract the first two colors\n",
+ "color1 = deep_palette[0]\n",
+ "color2 = deep_palette[1]\n",
+ "\n",
+ "#, notch=True\n",
+ "sns.set(style=\"white\", font_scale=1.5, font=\"sans-serif\", palette=\"Set1\")\n",
+ "\n",
+ "# get a label for outlier,figure out why there is error bar \n",
+ "# customize which counts as outlier\n",
+ "g = sns.catplot(x='value', y='balance', data=df, hue='treatment', kind='box',\n",
+ " hue_order=['reference', 'SLES'], n_boot = 1000,\n",
+ " height = 10, aspect = 1.8, legend = False, whis = [0,100])\n",
+ "g.set_xlabels('log fold change', fontsize=14)\n",
+ "g.set_ylabels('', fontsize=14)\n",
+ "plt.legend(fontsize = 20, frameon=False)\n",
+ "# plt.legend(loc='lower right', labels=['Control', 'Treatment'])\n",
+ "plt.legend(fontsize=20, frameon=False, loc='lower right')\n",
+ "plt.xticks(fontsize=12)\n",
+ "\n",
+ "ax = plt.gca()\n",
+ "leg = ax.legend(fontsize=20, frameon=False, loc='lower right')\n",
+ "leg.legendHandles[0].set_linewidth(8)\n",
+ "leg.legendHandles[0].set_color(color1)\n",
+ "leg.legendHandles[1].set_linewidth(8)\n",
+ "leg.legendHandles[1].set_color(color2)\n",
+ "leg.get_frame().set_linewidth(0)\n",
+ "#plt.xticks(rotation=60)\n",
+ "\n",
+ "# g.savefig('./out_sles/fig_sles.svg', format='svg')\n",
+ "stats.to_csv('./out_c_afermentans/stats_invitro_sles.csv')\n",
+ "stats"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "ef86565d-5578-45fd-8423-4ad797746748",
+ "metadata": {},
+ "source": [
+ "def paired_t_test(category, md, df):\n",
+ " \"\"\" Performs paired t-test on a dataframe with complementing\n",
+ " sample metadata.\n",
+ "\n",
+ " Parameters\n",
+ " ----------\n",
+ " category : str\n",
+ " Category to test in the dataframe df\n",
+ " md : pd.DataFrame\n",
+ " Sample metadata specifying groupings\n",
+ " df : pd.Dataframe\n",
+ " Counts of microbes across all samples\n",
+ " \n",
+ " Returns\n",
+ " -------\n",
+ " T : float\n",
+ " T-statistic\n",
+ " prob : float\n",
+ " p-value\n",
+ " lower : float\n",
+ " lower 95% CI bound\n",
+ " upper : float\n",
+ " upper 95% CI bound\n",
+ " ddof : int\n",
+ " degrees of freedom\n",
+ " \"\"\"\n",
+ " cats = list(md[category].value_counts().index)\n",
+ " src = md[md[category] == cats[0]].index\n",
+ " print(src)\n",
+ " dest = md[md[category] == cats[1]].index\n",
+ " return ttest_rel(df.loc[dest].values, df.loc[src].values)# RL"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 44,
+ "id": "ac689ee4-cdd1-4996-a62a-f96b95fefbbf",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def paired_t_test(category, md, df):\n",
+ " \"\"\" Performs paired t-test on a dataframe with complementing\n",
+ " sample metadata.\n",
+ "\n",
+ " Parameters\n",
+ " ----------\n",
+ " category : str\n",
+ " Category to test in the dataframe df\n",
+ " md : pd.DataFrame\n",
+ " Sample metadata specifying groupings\n",
+ " df : pd.Dataframe\n",
+ " Counts of microbes across all samples\n",
+ " \n",
+ " Returns\n",
+ " -------\n",
+ " T : float\n",
+ " T-statistic\n",
+ " prob : float\n",
+ " p-value\n",
+ " lower : float\n",
+ " lower 95% CI bound\n",
+ " upper : float\n",
+ " upper 95% CI bound\n",
+ " ddof : int\n",
+ " degrees of freedom\n",
+ " \"\"\"\n",
+ " cats = list(md[category].value_counts().index)\n",
+ " src = md[md[category] == cats[0]].index\n",
+ " print('***')\n",
+ " print(src)\n",
+ " dest = md[md[category] == cats[1]].index\n",
+ " return ttest_rel(df.loc[dest].values, df.loc[src].values)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 45,
+ "id": "43dcf082-3703-40f0-b5af-ca7175ace397",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "***\n",
+ "Index(['Comm-1_S4_L001', 'Comm-1_S4_L002', 'Comm-2_S14_L001',\n",
+ " 'Comm-2_S14_L002', 'Comm-3_S23_L001', 'Comm-3_S23_L002',\n",
+ " 'Comm-4_S32_L001', 'Comm-4_S32_L002'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['Comm-1_S4_L001', 'Comm-1_S4_L002', 'Comm-2_S14_L001',\n",
+ " 'Comm-2_S14_L002', 'Comm-3_S23_L001', 'Comm-3_S23_L002',\n",
+ " 'Comm-4_S32_L001', 'Comm-4_S32_L002'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['Comm-1_S4_L001', 'Comm-1_S4_L002', 'Comm-2_S14_L001',\n",
+ " 'Comm-2_S14_L002', 'Comm-3_S23_L001', 'Comm-3_S23_L002',\n",
+ " 'Comm-4_S32_L001', 'Comm-4_S32_L002'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['Comm-1_S4_L001', 'Comm-1_S4_L002', 'Comm-2_S14_L001',\n",
+ " 'Comm-2_S14_L002', 'Comm-3_S23_L001', 'Comm-3_S23_L002',\n",
+ " 'Comm-4_S32_L001', 'Comm-4_S32_L002'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['Comm-1_S4_L001', 'Comm-1_S4_L002', 'Comm-2_S14_L001',\n",
+ " 'Comm-2_S14_L002', 'Comm-3_S23_L001', 'Comm-3_S23_L002',\n",
+ " 'Comm-4_S32_L001', 'Comm-4_S32_L002'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['Comm-1_S4_L001', 'Comm-1_S4_L002', 'Comm-2_S14_L001',\n",
+ " 'Comm-2_S14_L002', 'Comm-3_S23_L001', 'Comm-3_S23_L002',\n",
+ " 'Comm-4_S32_L001', 'Comm-4_S32_L002'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['Comm-1_S4_L001', 'Comm-1_S4_L002', 'Comm-2_S14_L001',\n",
+ " 'Comm-2_S14_L002', 'Comm-3_S23_L001', 'Comm-3_S23_L002',\n",
+ " 'Comm-4_S32_L001', 'Comm-4_S32_L002'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['Comm-1_S4_L001', 'Comm-1_S4_L002', 'Comm-2_S14_L001',\n",
+ " 'Comm-2_S14_L002', 'Comm-3_S23_L001', 'Comm-3_S23_L002',\n",
+ " 'Comm-4_S32_L001', 'Comm-4_S32_L002'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['Comm-1_S4_L001', 'Comm-1_S4_L002', 'Comm-2_S14_L001',\n",
+ " 'Comm-2_S14_L002', 'Comm-3_S23_L001', 'Comm-3_S23_L002',\n",
+ " 'Comm-4_S32_L001', 'Comm-4_S32_L002'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['Comm-1_S4_L001', 'Comm-1_S4_L002', 'Comm-2_S14_L001',\n",
+ " 'Comm-2_S14_L002', 'Comm-3_S23_L001', 'Comm-3_S23_L002',\n",
+ " 'Comm-4_S32_L001', 'Comm-4_S32_L002'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['Comm-1_S4_L001', 'Comm-1_S4_L002', 'Comm-2_S14_L001',\n",
+ " 'Comm-2_S14_L002', 'Comm-3_S23_L001', 'Comm-3_S23_L002',\n",
+ " 'Comm-4_S32_L001', 'Comm-4_S32_L002'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['Comm-1_S4_L001', 'Comm-1_S4_L002', 'Comm-2_S14_L001',\n",
+ " 'Comm-2_S14_L002', 'Comm-3_S23_L001', 'Comm-3_S23_L002',\n",
+ " 'Comm-4_S32_L001', 'Comm-4_S32_L002'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['Comm-1_S4_L001', 'Comm-1_S4_L002', 'Comm-2_S14_L001',\n",
+ " 'Comm-2_S14_L002', 'Comm-3_S23_L001', 'Comm-3_S23_L002',\n",
+ " 'Comm-4_S32_L001', 'Comm-4_S32_L002'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['Comm-1_S4_L001', 'Comm-1_S4_L002', 'Comm-2_S14_L001',\n",
+ " 'Comm-2_S14_L002', 'Comm-3_S23_L001', 'Comm-3_S23_L002',\n",
+ " 'Comm-4_S32_L001', 'Comm-4_S32_L002'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['Comm-1_S4_L001', 'Comm-1_S4_L002', 'Comm-2_S14_L001',\n",
+ " 'Comm-2_S14_L002', 'Comm-3_S23_L001', 'Comm-3_S23_L002',\n",
+ " 'Comm-4_S32_L001', 'Comm-4_S32_L002'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['Comm-1_S4_L001', 'Comm-1_S4_L002', 'Comm-2_S14_L001',\n",
+ " 'Comm-2_S14_L002', 'Comm-3_S23_L001', 'Comm-3_S23_L002',\n",
+ " 'Comm-4_S32_L001', 'Comm-4_S32_L002'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['Comm-1_S4_L001', 'Comm-1_S4_L002', 'Comm-2_S14_L001',\n",
+ " 'Comm-2_S14_L002', 'Comm-3_S23_L001', 'Comm-3_S23_L002',\n",
+ " 'Comm-4_S32_L001', 'Comm-4_S32_L002'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['Comm-1_S4_L001', 'Comm-1_S4_L002', 'Comm-2_S14_L001',\n",
+ " 'Comm-2_S14_L002', 'Comm-3_S23_L001', 'Comm-3_S23_L002',\n",
+ " 'Comm-4_S32_L001', 'Comm-4_S32_L002'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['Comm-1_S4_L001', 'Comm-1_S4_L002', 'Comm-2_S14_L001',\n",
+ " 'Comm-2_S14_L002', 'Comm-3_S23_L001', 'Comm-3_S23_L002',\n",
+ " 'Comm-4_S32_L001', 'Comm-4_S32_L002'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['Comm-1_S4_L001', 'Comm-1_S4_L002', 'Comm-2_S14_L001',\n",
+ " 'Comm-2_S14_L002', 'Comm-3_S23_L001', 'Comm-3_S23_L002',\n",
+ " 'Comm-4_S32_L001', 'Comm-4_S32_L002'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['Comm-1_S4_L001', 'Comm-1_S4_L002', 'Comm-2_S14_L001',\n",
+ " 'Comm-2_S14_L002', 'Comm-3_S23_L001', 'Comm-3_S23_L002',\n",
+ " 'Comm-4_S32_L001', 'Comm-4_S32_L002'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['Comm-1_S4_L001', 'Comm-1_S4_L002', 'Comm-2_S14_L001',\n",
+ " 'Comm-2_S14_L002', 'Comm-3_S23_L001', 'Comm-3_S23_L002',\n",
+ " 'Comm-4_S32_L001', 'Comm-4_S32_L002'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['Comm-1_S4_L001', 'Comm-1_S4_L002', 'Comm-2_S14_L001',\n",
+ " 'Comm-2_S14_L002', 'Comm-3_S23_L001', 'Comm-3_S23_L002',\n",
+ " 'Comm-4_S32_L001', 'Comm-4_S32_L002'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['Comm-1_S4_L001', 'Comm-1_S4_L002', 'Comm-2_S14_L001',\n",
+ " 'Comm-2_S14_L002', 'Comm-3_S23_L001', 'Comm-3_S23_L002',\n",
+ " 'Comm-4_S32_L001', 'Comm-4_S32_L002'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['Comm-1_S4_L001', 'Comm-1_S4_L002', 'Comm-2_S14_L001',\n",
+ " 'Comm-2_S14_L002', 'Comm-3_S23_L001', 'Comm-3_S23_L002',\n",
+ " 'Comm-4_S32_L001', 'Comm-4_S32_L002'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['Comm-1_S4_L001', 'Comm-1_S4_L002', 'Comm-2_S14_L001',\n",
+ " 'Comm-2_S14_L002', 'Comm-3_S23_L001', 'Comm-3_S23_L002',\n",
+ " 'Comm-4_S32_L001', 'Comm-4_S32_L002'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['Comm-1_S4_L001', 'Comm-1_S4_L002', 'Comm-2_S14_L001',\n",
+ " 'Comm-2_S14_L002', 'Comm-3_S23_L001', 'Comm-3_S23_L002',\n",
+ " 'Comm-4_S32_L001', 'Comm-4_S32_L002'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['Comm-1_S4_L001', 'Comm-1_S4_L002', 'Comm-2_S14_L001',\n",
+ " 'Comm-2_S14_L002', 'Comm-3_S23_L001', 'Comm-3_S23_L002',\n",
+ " 'Comm-4_S32_L001', 'Comm-4_S32_L002'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['Comm-1_S4_L001', 'Comm-1_S4_L002', 'Comm-2_S14_L001',\n",
+ " 'Comm-2_S14_L002', 'Comm-3_S23_L001', 'Comm-3_S23_L002',\n",
+ " 'Comm-4_S32_L001', 'Comm-4_S32_L002'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['Comm-1_S4_L001', 'Comm-1_S4_L002', 'Comm-2_S14_L001',\n",
+ " 'Comm-2_S14_L002', 'Comm-3_S23_L001', 'Comm-3_S23_L002',\n",
+ " 'Comm-4_S32_L001', 'Comm-4_S32_L002'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['Comm-1_S4_L001', 'Comm-1_S4_L002', 'Comm-2_S14_L001',\n",
+ " 'Comm-2_S14_L002', 'Comm-3_S23_L001', 'Comm-3_S23_L002',\n",
+ " 'Comm-4_S32_L001', 'Comm-4_S32_L002'],\n",
+ " dtype='object', name='sampleid')\n",
+ "***\n",
+ "Index(['Comm-1_S4_L001', 'Comm-1_S4_L002', 'Comm-2_S14_L001',\n",
+ " 'Comm-2_S14_L002', 'Comm-3_S23_L001', 'Comm-3_S23_L002',\n",
+ " 'Comm-4_S32_L001', 'Comm-4_S32_L002'],\n",
+ " dtype='object', name='sampleid')\n"
+ ]
+ },
+ {
+ "data": {
+ "text/html": [
+ "\n",
+ "\n",
+ "
\n",
+ " \n",
+ " \n",
+ " | \n",
+ " 24hControl_vs_24hChemical_tstat | \n",
+ " 24hControl__vs_24hChemical_pvalue | \n",
+ " 24hControl_vs_24hChemical_upperCI | \n",
+ " 24hControl_vs_24hChemical_lowerCI | \n",
+ "
\n",
+ " \n",
+ " \n",
+ " \n",
+ " | log(C.acnes/C.afermentans) | \n",
+ " -5.159104 | \n",
+ " 0.001311 | \n",
+ " -4.112857 | \n",
+ " -1.527606 | \n",
+ "
\n",
+ " \n",
+ " | log(M.luteus/C.afermentans) | \n",
+ " 0.827653 | \n",
+ " 0.435172 | \n",
+ " -0.027550 | \n",
+ " 0.057221 | \n",
+ "
\n",
+ " \n",
+ " | log(S.aureus/C.afermentans) | \n",
+ " -9.432206 | \n",
+ " 0.000031 | \n",
+ " -4.040212 | \n",
+ " -2.420525 | \n",
+ "
\n",
+ " \n",
+ " | log(S.capitis/C.afermentans) | \n",
+ " -3.635776 | \n",
+ " 0.008334 | \n",
+ " -0.396925 | \n",
+ " -0.084086 | \n",
+ "
\n",
+ " \n",
+ " | log(S.epidermidis/C.afermentans) | \n",
+ " 2.843834 | \n",
+ " 0.024907 | \n",
+ " 0.082282 | \n",
+ " 0.894311 | \n",
+ "
\n",
+ " \n",
+ " | log(S.hominis/C.afermentans) | \n",
+ " 3.933570 | \n",
+ " 0.005650 | \n",
+ " 0.197734 | \n",
+ " 0.793762 | \n",
+ "
\n",
+ " \n",
+ " | log(S.warneri/C.afermentans) | \n",
+ " 2.485341 | \n",
+ " 0.041881 | \n",
+ " 0.017821 | \n",
+ " 0.715997 | \n",
+ "
\n",
+ " \n",
+ " | log(S.mitis/C.afermentans) | \n",
+ " 1.308089 | \n",
+ " 0.232158 | \n",
+ " -0.234270 | \n",
+ " 0.814365 | \n",
+ "
\n",
+ " \n",
+ "
\n",
+ "