Code13_DeltaBMI-misclassification_gut-microbiome.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Multiomics BMI Paper — Relationships between ∆BMI-based Misclassification and Gut Microbiome Alpha-diversity"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***Analyzed by Tomasz Wilmanski originally, and modified by Kengo Watanabe***  \n",
    "\n",
    "This Jupyter Notebook (with Python 3 kernel) assessed the relationships between ∆BMI-based misclassification (i.e., BMI class vs. biological BMI class) and gut microbiome alpha-diversity metrics using regression analysis (in the baseline Arivale cohort).  \n",
    "\n",
    "Input files:  \n",
    "* Arivale baseline biological BMIs and covariates: 220803_Multiomics-BMI-NatMed1stRevision_DeltaBMI-misclassification_biologicalBMI-baseline-summary-BothSex.tsv  \n",
    "* Arivale baseline gut microbiome alpha-diversity metrics: 220902_Multiomics-BMI-NatMed1stRevision_Microbiome-DataCleaning_AlphaDiversity-and-TaxonAbundance_final.tsv  \n",
    "\n",
    "Output figures and tables:  \n",
    "* Figure 4b  \n",
    "* Table for Supplementary Data 6  \n",
    "\n",
    "Original notebook (memo for my future tracing):  \n",
    "* wenceslaus:\\[JupyterLab HOME\\]/220621_Multiomics-BMI-NatMedRevision/220905_Multiomics-BMI-NatMed1stRevision_Microbiome-DeltaBMI-wenceslaus.ipynb  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "# packages in environment at /opt/conda/envs/arivale-py3:\n",
      "#\n",
      "# Name                    Version                   Build  Channel\n",
      "_libgcc_mutex             0.1                 conda_forge    conda-forge\n",
      "_openmp_mutex             4.5                       1_gnu    conda-forge\n",
      "analytics                 0.1                      pypi_0    pypi\n",
      "argon2-cffi               21.1.0           py39h3811e60_0    conda-forge\n",
      "arivale-data-interface    0.1.0                    pypi_0    pypi\n",
      "async_generator           1.10                       py_0    conda-forge\n",
      "atk-1.0                   2.36.0               h3371d22_4    conda-forge\n",
      "attrs                     21.2.0             pyhd8ed1ab_0    conda-forge\n",
      "backcall                  0.2.0              pyh9f0ad1d_0    conda-forge\n",
      "backports                 1.0                        py_2    conda-forge\n",
      "backports.functools_lru_cache 1.6.4              pyhd8ed1ab_0    conda-forge\n",
      "biopython                 1.79             py39h3811e60_0    conda-forge\n",
      "bleach                    4.1.0              pyhd8ed1ab_0    conda-forge\n",
      "bokeh                     2.4.1            py39hf3d152e_1    conda-forge\n",
      "boto                      2.49.0                     py_0    conda-forge\n",
      "boto3                     1.19.2                   pypi_0    pypi\n",
      "botocore                  1.22.2                   pypi_0    pypi\n",
      "brotlipy                  0.7.0           py39h3811e60_1001    conda-forge\n",
      "bzip2                     1.0.8                h7f98852_4    conda-forge\n",
      "c-ares                    1.17.2               h7f98852_0    conda-forge\n",
      "ca-certificates           2022.6.15.1          ha878542_0    conda-forge\n",
      "cachecontrol              0.12.6                     py_0    conda-forge\n",
      "cairo                     1.16.0            h6cf1ce9_1008    conda-forge\n",
      "certifi                   2022.6.15.1        pyhd8ed1ab_0    conda-forge\n",
      "cffi                      1.14.6           py39h4bc2ebd_1    conda-forge\n",
      "chardet                   4.0.0            py39hf3d152e_1    conda-forge\n",
      "charset-normalizer        2.0.0              pyhd8ed1ab_0    conda-forge\n",
      "click                     8.0.3            py39hf3d152e_0    conda-forge\n",
      "colorama                  0.4.4              pyh9f0ad1d_0    conda-forge\n",
      "cryptography              35.0.0           py39h95dcef6_1    conda-forge\n",
      "cycler                    0.10.0                     py_2    conda-forge\n",
      "cython                    0.29.24          py39he80948d_0    conda-forge\n",
      "debugpy                   1.4.1            py39he80948d_0    conda-forge\n",
      "decorator                 4.4.2                      py_0    conda-forge\n",
      "defusedxml                0.7.1              pyhd8ed1ab_0    conda-forge\n",
      "english                   2020.7.0           pyhd8ed1ab_0    conda-forge\n",
      "entrypoints               0.3             py39hde42818_1002    conda-forge\n",
      "et_xmlfile                1.0.1                   py_1001    conda-forge\n",
      "expat                     2.4.1                h9c3ff4c_0    conda-forge\n",
      "font-ttf-dejavu-sans-mono 2.37                 hab24e00_0    conda-forge\n",
      "font-ttf-inconsolata      3.000                h77eed37_0    conda-forge\n",
      "font-ttf-source-code-pro  2.038                h77eed37_0    conda-forge\n",
      "font-ttf-ubuntu           0.83                 hab24e00_0    conda-forge\n",
      "fontconfig                2.13.1            hba837de_1005    conda-forge\n",
      "fonts-conda-ecosystem     1                             0    conda-forge\n",
      "fonts-conda-forge         1                             0    conda-forge\n",
      "freetype                  2.10.4               h0708190_1    conda-forge\n",
      "fribidi                   1.0.10               h36c2ea0_0    conda-forge\n",
      "future                    0.18.2           py39hf3d152e_3    conda-forge\n",
      "gdk-pixbuf                2.42.6               h04a7f16_0    conda-forge\n",
      "gettext                   0.19.8.1          h73d1719_1008    conda-forge\n",
      "giflib                    5.2.1                h36c2ea0_2    conda-forge\n",
      "graphite2                 1.3.13            h58526e2_1001    conda-forge\n",
      "graphviz                  2.49.1               h85b4f2f_0    conda-forge\n",
      "gtk2                      2.24.33              h539f30e_1    conda-forge\n",
      "gts                       0.7.6                h64030ff_2    conda-forge\n",
      "harfbuzz                  3.0.0                h83ec7ef_1    conda-forge\n",
      "hdmedians                 0.14.2           py39hce5d2b2_0    conda-forge\n",
      "icu                       68.1                 h58526e2_0    conda-forge\n",
      "idna                      3.1                pyhd3deb0d_0    conda-forge\n",
      "importlib-metadata        4.8.1            py39hf3d152e_0    conda-forge\n",
      "iniconfig                 1.1.1              pyh9f0ad1d_0    conda-forge\n",
      "interchange               2021.0.4           pyhd8ed1ab_0    conda-forge\n",
      "ipykernel                 6.4.2            py39hef51801_0    conda-forge\n",
      "ipython                   7.28.0           py39hef51801_0    conda-forge\n",
      "ipython_genutils          0.2.0                      py_1    conda-forge\n",
      "ipywidgets                7.6.5              pyhd8ed1ab_0    conda-forge\n",
      "jbig                      2.1               h7f98852_2003    conda-forge\n",
      "jedi                      0.18.0           py39hf3d152e_2    conda-forge\n",
      "jinja2                    3.0.2              pyhd8ed1ab_0    conda-forge\n",
      "jmespath                  0.10.0                   pypi_0    pypi\n",
      "joblib                    1.1.0              pyhd8ed1ab_0    conda-forge\n",
      "jpeg                      9d                   h36c2ea0_0    conda-forge\n",
      "jsonschema                4.1.2              pyhd8ed1ab_0    conda-forge\n",
      "jupyter_client            7.0.6              pyhd8ed1ab_0    conda-forge\n",
      "jupyter_core              4.8.1            py39hf3d152e_0    conda-forge\n",
      "jupyterlab_pygments       0.1.2              pyh9f0ad1d_0    conda-forge\n",
      "jupyterlab_widgets        1.0.2              pyhd8ed1ab_0    conda-forge\n",
      "kiwisolver                1.3.2            py39h1a9c180_0    conda-forge\n",
      "krb5                      1.19.2               hcc1bbae_2    conda-forge\n",
      "lcms2                     2.12                 hddcbb42_0    conda-forge\n",
      "ld_impl_linux-64          2.36.1               hea4e1c9_2    conda-forge\n",
      "lerc                      2.2.1                h9c3ff4c_0    conda-forge\n",
      "libblas                   3.9.0           12_linux64_openblas    conda-forge\n",
      "libcblas                  3.9.0           12_linux64_openblas    conda-forge\n",
      "libcurl                   7.79.1               h2574ce0_1    conda-forge\n",
      "libdeflate                1.7                  h7f98852_5    conda-forge\n",
      "libedit                   3.1.20191231         he28a2e2_2    conda-forge\n",
      "libev                     4.33                 h516909a_1    conda-forge\n",
      "libffi                    3.4.2                h9c3ff4c_4    conda-forge\n",
      "libgcc-ng                 11.2.0              h1d223b6_11    conda-forge\n",
      "libgd                     2.3.3                h6ad9fb6_0    conda-forge\n",
      "libgfortran-ng            11.2.0              h69a702a_11    conda-forge\n",
      "libgfortran5              11.2.0              h5c6108e_11    conda-forge\n",
      "libglib                   2.70.0               h174f98d_1    conda-forge\n",
      "libgomp                   11.2.0              h1d223b6_11    conda-forge\n",
      "libiconv                  1.16                 h516909a_0    conda-forge\n",
      "liblapack                 3.9.0           12_linux64_openblas    conda-forge\n",
      "libnghttp2                1.43.0               h812cca2_1    conda-forge\n",
      "libopenblas               0.3.18          pthreads_h8fe5266_0    conda-forge\n",
      "libpng                    1.6.37               h21135ba_2    conda-forge\n",
      "libpq                     13.3                 hd57d9b9_1    conda-forge\n",
      "librsvg                   2.52.2               hc3c00ef_0    conda-forge\n",
      "libsodium                 1.0.18               h36c2ea0_1    conda-forge\n",
      "libssh2                   1.10.0               ha56f1ee_2    conda-forge\n",
      "libstdcxx-ng              11.2.0              he4da1e4_11    conda-forge\n",
      "libtiff                   4.3.0                hf544144_1    conda-forge\n",
      "libtool                   2.4.6             h9c3ff4c_1008    conda-forge\n",
      "libuuid                   2.32.1            h7f98852_1000    conda-forge\n",
      "libwebp                   1.2.1                h3452ae3_0    conda-forge\n",
      "libwebp-base              1.2.1                h7f98852_0    conda-forge\n",
      "libxcb                    1.13              h7f98852_1003    conda-forge\n",
      "libxml2                   2.9.12               h72842e0_0    conda-forge\n",
      "libxslt                   1.1.33               h15afd5d_2    conda-forge\n",
      "libzlib                   1.2.11            h36c2ea0_1013    conda-forge\n",
      "lockfile                  0.12.2                     py_1    conda-forge\n",
      "lxml                      4.6.3            py39h107f48f_0    conda-forge\n",
      "lz4-c                     1.9.3                h9c3ff4c_1    conda-forge\n",
      "markupsafe                2.0.1            py39h3811e60_0    conda-forge\n",
      "matplotlib-base           3.4.3            py39h2fa2bec_1    conda-forge\n",
      "matplotlib-inline         0.1.3              pyhd8ed1ab_0    conda-forge\n",
      "missingpy                 0.2.0                    pypi_0    pypi\n",
      "mistune                   0.8.4           py39h3811e60_1004    conda-forge\n",
      "more-itertools            8.10.0             pyhd8ed1ab_0    conda-forge\n",
      "mscorefonts               0.0.1                         3    conda-forge\n",
      "msgpack-python            1.0.2            py39h1a9c180_1    conda-forge\n",
      "natsort                   7.1.1              pyhd8ed1ab_0    conda-forge\n",
      "nbclient                  0.5.4              pyhd8ed1ab_0    conda-forge\n",
      "nbconvert                 6.2.0            py39hf3d152e_0    conda-forge\n",
      "nbformat                  5.1.3              pyhd8ed1ab_0    conda-forge\n",
      "ncurses                   6.2                  h58526e2_4    conda-forge\n",
      "neo4j-python-driver       4.3.7              pyhd8ed1ab_0    conda-forge\n",
      "neobolt                   1.7.17           py39h3811e60_2    conda-forge\n",
      "neotime                   1.7.4                      py_0    conda-forge\n",
      "nest-asyncio              1.5.1              pyhd8ed1ab_0    conda-forge\n",
      "networkx                  2.6.3              pyhd8ed1ab_0    conda-forge\n",
      "notebook                  6.4.5              pyha770c72_0    conda-forge\n",
      "numpy                     1.21.3           py39hdbf815f_0    conda-forge\n",
      "olefile                   0.46               pyh9f0ad1d_1    conda-forge\n",
      "openjpeg                  2.4.0                hb52868f_1    conda-forge\n",
      "openpyxl                  3.0.9              pyhd8ed1ab_0    conda-forge\n",
      "openssl                   1.1.1l               h7f98852_0    conda-forge\n",
      "packaging                 21.0               pyhd8ed1ab_0    conda-forge\n",
      "pandas                    1.3.4            py39hde0f152_0    conda-forge\n",
      "pandoc                    2.15                 h7f98852_0    conda-forge\n",
      "pandocfilters             1.5.0              pyhd8ed1ab_0    conda-forge\n",
      "pango                     1.48.10              h54213e6_2    conda-forge\n",
      "parso                     0.8.2              pyhd8ed1ab_0    conda-forge\n",
      "patsy                     0.5.2              pyhd8ed1ab_0    conda-forge\n",
      "pcre                      8.45                 h9c3ff4c_0    conda-forge\n",
      "pexpect                   4.8.0              pyh9f0ad1d_2    conda-forge\n",
      "pickleshare               0.7.5           py39hde42818_1002    conda-forge\n",
      "pillow                    8.3.2            py39ha612740_0    conda-forge\n",
      "pip                       21.3.1             pyhd8ed1ab_0    conda-forge\n",
      "pixman                    0.40.0               h36c2ea0_0    conda-forge\n",
      "plotly                    5.3.1              pyhd8ed1ab_0    conda-forge\n",
      "pluggy                    1.0.0            py39hf3d152e_1    conda-forge\n",
      "prometheus_client         0.11.0             pyhd8ed1ab_0    conda-forge\n",
      "prompt-toolkit            3.0.21             pyha770c72_0    conda-forge\n",
      "prompt_toolkit            3.0.21               hd8ed1ab_0    conda-forge\n",
      "psycopg2                  2.9.1            py39h3811e60_0    conda-forge\n",
      "pthread-stubs             0.4               h36c2ea0_1001    conda-forge\n",
      "ptyprocess                0.7.0              pyhd3deb0d_0    conda-forge\n",
      "py                        1.10.0             pyhd3deb0d_0    conda-forge\n",
      "py2neo                    2021.2.3           pyhd8ed1ab_0    conda-forge\n",
      "pycparser                 2.20               pyh9f0ad1d_2    conda-forge\n",
      "pygments                  2.10.0             pyhd8ed1ab_0    conda-forge\n",
      "pygraphviz                1.7              py39h78163bd_0    conda-forge\n",
      "pyopenssl                 21.0.0             pyhd8ed1ab_0    conda-forge\n",
      "pyparsing                 2.4.7              pyh9f0ad1d_0    conda-forge\n",
      "pyrsistent                0.17.3           py39h3811e60_2    conda-forge\n",
      "pysam                     0.17.0           py39h051187c_0    bioconda\n",
      "pysocks                   1.7.1            py39hf3d152e_3    conda-forge\n",
      "pytabix                   0.1              py39h98c8e45_1    bioconda\n",
      "pytest                    6.2.5            py39hf3d152e_0    conda-forge\n",
      "python                    3.9.7           hb7a2778_3_cpython    conda-forge\n",
      "python-dateutil           2.8.2              pyhd8ed1ab_0    conda-forge\n",
      "python-graphviz           0.17               pyhaef67bd_0    conda-forge\n",
      "python_abi                3.9                      2_cp39    conda-forge\n",
      "pytz                      2021.3             pyhd8ed1ab_0    conda-forge\n",
      "pyvcf                     0.6.8           py39hde42818_1002    conda-forge\n",
      "pyyaml                    6.0              py39h3811e60_0    conda-forge\n",
      "pyzmq                     22.3.0           py39h37b5a0c_0    conda-forge\n",
      "readline                  8.1                  h46c0cb4_0    conda-forge\n",
      "requests                  2.26.0             pyhd8ed1ab_0    conda-forge\n",
      "s3transfer                0.5.0                    pypi_0    pypi\n",
      "scikit-bio                0.5.6            py39h16ac069_4    conda-forge\n",
      "scikit-learn              1.0.1            py39h7c5d8c9_0    conda-forge\n",
      "scipy                     1.7.1            py39hee8e79c_0    conda-forge\n",
      "seaborn                   0.11.2               hd8ed1ab_0    conda-forge\n",
      "seaborn-base              0.11.2             pyhd8ed1ab_0    conda-forge\n",
      "send2trash                1.8.0              pyhd8ed1ab_0    conda-forge\n",
      "setuptools                58.2.0           py39hf3d152e_0    conda-forge\n",
      "six                       1.16.0             pyh6c4a22f_0    conda-forge\n",
      "sqlite                    3.36.0               h9cd32fc_2    conda-forge\n",
      "statsmodels               0.13.0           py39hce5d2b2_0    conda-forge\n",
      "tenacity                  8.0.1              pyhd8ed1ab_0    conda-forge\n",
      "terminado                 0.12.1           py39hf3d152e_0    conda-forge\n",
      "testpath                  0.5.0              pyhd8ed1ab_0    conda-forge\n",
      "threadpoolctl             3.0.0              pyh8a188c0_0    conda-forge\n",
      "tk                        8.6.11               h27826a3_1    conda-forge\n",
      "toml                      0.10.2             pyhd8ed1ab_0    conda-forge\n",
      "tornado                   6.1              py39h3811e60_1    conda-forge\n",
      "traitlets                 5.1.0              pyhd8ed1ab_0    conda-forge\n",
      "typing_extensions         3.10.0.2           pyha770c72_0    conda-forge\n",
      "tzdata                    2021e                he74cb21_0    conda-forge\n",
      "urllib3                   1.26.7             pyhd8ed1ab_0    conda-forge\n",
      "wcwidth                   0.2.5              pyh9f0ad1d_2    conda-forge\n",
      "webencodings              0.5.1                      py_1    conda-forge\n",
      "wheel                     0.37.0             pyhd8ed1ab_1    conda-forge\n",
      "widgetsnbextension        3.5.1            py39hf3d152e_4    conda-forge\n",
      "wordcloud                 1.8.1                    pypi_0    pypi\n",
      "xlrd                      2.0.1              pyhd8ed1ab_3    conda-forge\n",
      "xorg-kbproto              1.0.7             h7f98852_1002    conda-forge\n",
      "xorg-libice               1.0.10               h7f98852_0    conda-forge\n",
      "xorg-libsm                1.2.3             hd9c2040_1000    conda-forge\n",
      "xorg-libx11               1.7.2                h7f98852_0    conda-forge\n",
      "xorg-libxau               1.0.9                h7f98852_0    conda-forge\n",
      "xorg-libxdmcp             1.1.3                h7f98852_0    conda-forge\n",
      "xorg-libxext              1.3.4                h7f98852_1    conda-forge\n",
      "xorg-libxrender           0.9.10            h7f98852_1003    conda-forge\n",
      "xorg-renderproto          0.11.1            h7f98852_1002    conda-forge\n",
      "xorg-xextproto            7.3.0             h7f98852_1002    conda-forge\n",
      "xorg-xproto               7.0.31            h7f98852_1007    conda-forge\n",
      "xz                        5.2.5                h516909a_1    conda-forge\n",
      "yaml                      0.2.5                h516909a_0    conda-forge\n",
      "zeromq                    4.3.4                h9c3ff4c_1    conda-forge\n",
      "zipp                      3.6.0              pyhd8ed1ab_0    conda-forge\n",
      "zlib                      1.2.11            h36c2ea0_1013    conda-forge\n",
      "zstd                      1.5.0                ha95c52a_0    conda-forge\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import scipy.stats as stats\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "import seaborn as sns\n",
    "#For Arial font\n",
    "#!conda install -c conda-forge -y mscorefonts\n",
    "##-> The below was also needed in matplotlib 3.4.2\n",
    "#import shutil\n",
    "#import matplotlib\n",
    "#shutil.rmtree(matplotlib.get_cachedir())\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "from IPython.display import display\n",
    "import time\n",
    "\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "import statsmodels.formula.api as smf\n",
    "from statsmodels.stats import multitest as multi\n",
    "from decimal import Decimal, ROUND_HALF_UP\n",
    "\n",
    "!conda list"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Prepare datasets"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> The necessary files were copied from the dalek server in advance.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1-1. Alpha-diversity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Import cleaned table for baseline gut microbiome data\n",
    "fileDir = './ImportData/'\n",
    "ipynbName = '220902_Multiomics-BMI-NatMed1stRevision_Microbiome-DataCleaning_'\n",
    "fileName = 'AlphaDiversity-and-TaxonAbundance_final.tsv'\n",
    "tempDF = pd.read_csv(fileDir+ipynbName+fileName, sep='\\t', dtype={'public_client_id': str})\n",
    "tempDF = tempDF.set_index('public_client_id')\n",
    "\n",
    "#Take only the alpha-diversity metrics\n",
    "tempL = ['Observed', 'Shannon', 'Chao1']\n",
    "tempDF = tempDF[tempL]\n",
    "\n",
    "display(tempDF)\n",
    "tempDF1 = tempDF.describe()\n",
    "tempDF1.loc['Skewness'] = stats.skew(tempDF)\n",
    "display(tempDF1)\n",
    "\n",
    "divDF = tempDF"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> –> Based on the Tom's original analysis, square transformation is applied to Shannon, while square root transformation is applied to Observed and Chao1.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempDF = divDF\n",
    "\n",
    "#Transform to less skewed distribution\n",
    "tempDF['Observed_tf'] = np.sqrt(tempDF['Observed'])\n",
    "tempDF['Shannon_tf'] = tempDF['Shannon']**2\n",
    "tempDF['Chao1_tf'] = np.sqrt(tempDF['Chao1'])\n",
    "\n",
    "display(tempDF)\n",
    "tempDF1 = tempDF.describe()\n",
    "tempDF1.loc['Skewness'] = stats.skew(tempDF)\n",
    "display(tempDF1)\n",
    "\n",
    "#Update\n",
    "divDF = tempDF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Check distribution and probability plot\n",
    "tempDF = divDF\n",
    "tempD1 = {'Observed':'tab:blue', 'Shannon':'tab:orange', 'Chao1':'tab:green',\n",
    "          'Observed_tf':'tab:blue', 'Shannon_tf':'tab:orange', 'Chao1_tf':'tab:green'}\n",
    "tempD2 = {'Observed':'Raw value', 'Shannon':'Raw value', 'Chao1':'Raw value',\n",
    "          'Observed_tf':'Transformed value', 'Shannon_tf':'Transformed value', 'Chao1_tf':'Transformed value'}\n",
    "\n",
    "sns.set(style='ticks', font='Arial', context='notebook')\n",
    "fig, axes = plt.subplots(nrows=2, ncols=6, figsize=(15, 5), sharex=False, sharey=False)\n",
    "for ax_i, ax in enumerate(axes.flat):\n",
    "    col_i = (ax_i//12)*6 + ax_i%6\n",
    "    if col_i<len(tempDF.columns):\n",
    "        col_n = tempDF.columns.tolist()[col_i]\n",
    "        if (ax_i//6)%2 == 0:\n",
    "            sns.set(style='ticks', font='Arial', context='notebook')\n",
    "            sns.distplot(tempDF[col_n], color=tempD1[col_n], ax=ax)\n",
    "            sns.despine()\n",
    "            ax.set_title(col_n, fontsize='large')\n",
    "            ax.set_xlabel(tempD2[col_n])\n",
    "        else:\n",
    "            sns.set(style='whitegrid', font='Arial', context='notebook')\n",
    "            stats.probplot(tempDF[col_n], plot=ax)\n",
    "            ax.get_lines()[0].set_markerfacecolor(tempD1[col_n])\n",
    "            ax.get_lines()[0].set_markeredgecolor(tempD1[col_n])\n",
    "            ax.get_lines()[1].set_color('k')\n",
    "            ax.get_lines()[1].set_linewidth(3)\n",
    "            skewness = stats.skew(tempDF[col_n])\n",
    "            ax.set_title('Skewness: '+str(round(skewness, 3)), fontsize='large')\n",
    "    else:\n",
    "        ax.plot([0, 1], [0, 1], color='k')\n",
    "        ax.plot([0, 1], [1, 0], color='k')\n",
    "fig.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1-2. ∆BMI-derived misclassification and covariates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Import cleaned table for baseline measured and biological BMIs\n",
    "fileDir = './ImportData/'\n",
    "ipynbName = '220803_Multiomics-BMI-NatMed1stRevision_DeltaBMI-misclassification_'\n",
    "fileName = 'biologicalBMI-baseline-summary-BothSex.tsv'\n",
    "tempDF = pd.read_csv(fileDir+ipynbName+fileName, sep='\\t', dtype={'public_client_id': str})\n",
    "tempDF = tempDF.set_index('public_client_id')\n",
    "print('Original:', len(tempDF))\n",
    "\n",
    "#Take the participants having alpha-diversity\n",
    "tempDF = tempDF.loc[tempDF.index.isin(divDF.index.tolist())]\n",
    "print(' -> with alpha-diversity:', len(tempDF))\n",
    "\n",
    "#Clean to handle easier in this notebook\n",
    "tempDF.columns = tempDF.columns.str.replace('Base', '')\n",
    "\n",
    "#Misclassification\n",
    "tempL = ['MetBMI', 'ProtBMI', 'ChemBMI', 'CombiBMI']\n",
    "for bbmi in tempL:\n",
    "    tempL1 = []\n",
    "    for row_i in range(len(tempDF)):\n",
    "        bmi_class = tempDF['BMI_class'].iloc[row_i]\n",
    "        bbmi_class = tempDF[bbmi+'_class'].iloc[row_i]\n",
    "        if bmi_class==bbmi_class:\n",
    "            tempL1.append('Matched')\n",
    "        else:\n",
    "            tempL1.append('Mismatched')\n",
    "    tempDF['vs_'+bbmi+'_class'] = tempL1\n",
    "\n",
    "display(tempDF)\n",
    "print('NaN in DF:', tempDF.isnull().to_numpy().sum(axis=None))\n",
    "\n",
    "bmiDF = tempDF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempDF = bmiDF\n",
    "\n",
    "#Check\n",
    "for bmi_class in ['Normal', 'Obese']:\n",
    "    print(bmi_class)\n",
    "    tempDF1 = tempDF.loc[tempDF['BMI_class']==bmi_class]\n",
    "    for bbmi in tempL:\n",
    "        display(tempDF1.groupby('vs_'+bbmi+'_class')['BMI'].describe())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> –> It would be safer to adjust the baseline BMI in the statistical tests.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempDF = bmiDF\n",
    "\n",
    "#Select the misclassification and covariates (just for the display in Jupyter notebook)\n",
    "tempL1 = tempDF.loc[:, tempDF.columns.str.contains('_class')].columns.tolist()\n",
    "tempL2 = ['BMI', 'Sex', 'Age', 'PC1', 'PC2', 'PC3', 'PC4', 'PC5']\n",
    "tempL = [col_n for sublist in [tempL1, tempL2] for col_n in sublist]\n",
    "tempDF1 = tempDF[tempL]\n",
    "\n",
    "display(tempDF1)\n",
    "\n",
    "#Update\n",
    "bmiDF = tempDF1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Regression analysis for the alpha-diversity metrics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Because the datasets are almost ready for use, only the simple processing steps are required for OLS linear regression:  \n",
    "> * ~~Missingness: To maximize the sample size for each regression, dropping NaN is performed after selecting metric.~~  \n",
    "> * Centering: standardization is applied to both dependent and independent variables (including covariates).  \n",
    ">\n",
    "> Hence, the remaining processing steps are implemented during the for-loop for each regression.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2-1. Perform OLS linear regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Model: Metric ~ b0 + b1\\*C(Misclassification) + b2\\*BMI + b3\\*C(Sex) + b4\\*Age + b5\\*AncestryPCs  \n",
    "> Main aim: Assess the difference in each metric between the matched and mismatched BMI class.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempDF1 = divDF.loc[:, divDF.columns.str.contains('_tf')]#Transformed metric values\n",
    "tempDF2 = bmiDF\n",
    "tempL1 = ['Normal', 'Obese']\n",
    "tempL2 = ['MetBMI', 'ProtBMI', 'ChemBMI', 'CombiBMI']\n",
    "\n",
    "t_start = time.time()\n",
    "tempD1 = {}\n",
    "for metric in tempDF1.columns.tolist():\n",
    "    tempD2 = {}\n",
    "    for bmi_class in tempL1:\n",
    "        #Processing for OLS linear regression\n",
    "        ##Gather all necessary variables into a single DF\n",
    "        tempS = tempDF1[metric]\n",
    "        tempDF = pd.merge(tempS, tempDF2, left_index=True, right_index=True, how='left')\n",
    "        ##Select the target participants\n",
    "        tempDF = tempDF.loc[tempDF['BMI_class']==bmi_class]\n",
    "        ##Drop NaN in the metric values\n",
    "        #tempDF = tempDF.dropna()\n",
    "        ##Z-score transformation\n",
    "        tempDF4 = tempDF.select_dtypes(include=[np.number])\n",
    "        scaler = StandardScaler(copy=True, with_mean=True, with_std=True)\n",
    "        tempA = scaler.fit_transform(tempDF4)#Column direction\n",
    "        tempDF4 = pd.DataFrame(data=tempA, index=tempDF4.index, columns=tempDF4.columns)\n",
    "        ###Recover the categorical variables\n",
    "        tempDF5 = tempDF.select_dtypes(exclude=[np.number])\n",
    "        tempDF = pd.merge(tempDF4, tempDF5, left_index=True, right_index=True, how='left')\n",
    "        ##Rename the dependent variable\n",
    "        tempDF = tempDF.rename(columns={metric:'Metric'})\n",
    "        ##Add a constant for the intercept\n",
    "        ###–> In statsmodels, a constant is automatically added as well as R!\n",
    "        \n",
    "        tempD3 = {}\n",
    "        for bbmi in tempL2:\n",
    "            #Processing specific to each bBMI\n",
    "            ##Sort to make bcoef = 0 and 1 for Matched and Mismatched\n",
    "            tempDF = tempDF.sort_values(by='vs_'+bbmi+'_class', ascending=True)\n",
    "            ##One-hot encoding for categorical covariates\n",
    "            ###–> In statsmodels, categorical variables are automatically recognized!\n",
    "            \n",
    "            #OLS linear regression\n",
    "            ##Fit univariate model\n",
    "            formula = 'Metric ~ C(vs_'+bbmi+'_class)'\n",
    "            fit_res1 = smf.ols(formula, data=tempDF).fit()\n",
    "            ##Fit full model\n",
    "            formula = 'Metric ~ C(vs_'+bbmi+'_class)'\\\n",
    "                '+ BMI + C(Sex) + Age + PC1 + PC2 + PC3 + PC4 + PC5'\n",
    "            fit_res2 = smf.ols(formula, data=tempDF).fit()\n",
    "            \n",
    "            #Summarize the result\n",
    "            tempS = pd.Series().astype('float64')\n",
    "            ##Save the sample size for each group\n",
    "            size1 = len(tempDF.loc[tempDF['vs_'+bbmi+'_class']=='Matched'])\n",
    "            size2 = len(tempDF.loc[tempDF['vs_'+bbmi+'_class']=='Mismatched'])\n",
    "            tempS.loc['vs'+bbmi+'class_nMatched'] = size1\n",
    "            tempS.loc['vs'+bbmi+'class_nMismatched'] = size2\n",
    "            ##Save R2 [%]\n",
    "            tempS.loc['vs'+bbmi+'class_UnivarR2'] = fit_res1.rsquared*100\n",
    "            tempS.loc['vs'+bbmi+'class_R2'] = fit_res2.rsquared*100\n",
    "            ##Save beta-coefficient of the target variable\n",
    "            tempS.loc['vs'+bbmi+'class_Bcoef'] = fit_res2.params['C(vs_'+bbmi+'_class)[T.Mismatched]']\n",
    "            tempS.loc['vs'+bbmi+'class_BcoefSE'] = fit_res2.bse['C(vs_'+bbmi+'_class)[T.Mismatched]']\n",
    "            ##Save t-statistic of the target variable\n",
    "            tempS.loc['vs'+bbmi+'class_tStat'] = fit_res2.tvalues['C(vs_'+bbmi+'_class)[T.Mismatched]']\n",
    "            ##Save P-value of the target variable\n",
    "            tempS.loc['vs'+bbmi+'class_Pval'] = fit_res2.pvalues['C(vs_'+bbmi+'_class)[T.Mismatched]']\n",
    "            ##Add dummy adjusted P-value rows for now\n",
    "            tempS.loc['vs'+bbmi+'class_AdjPval_within1'] = 1.0\n",
    "            tempS.loc['vs'+bbmi+'class_AdjPval_within2'] = 1.0\n",
    "            tempS.loc['vs'+bbmi+'class_AdjPval_all'] = 1.0\n",
    "            \n",
    "            tempD3[bbmi] = tempS\n",
    "        \n",
    "        #Clean the results (pd.Series) across bBMIs\n",
    "        ##Prepare common summary metrics: sample size, residual degrees of freedom\n",
    "        tempS1 = pd.Series().astype('float64')\n",
    "        tempS1.loc['N'] = len(tempDF)\n",
    "        tempS1.loc['DoF'] = int(fit_res2.df_resid)#Use the last result object but same b/w bBMIs\n",
    "        ##Combine each result\n",
    "        tempS2 = pd.concat(list(tempD3.values()), axis=0)\n",
    "        tempS = pd.concat([tempS1, tempS2], axis=0)\n",
    "        ##Convert to DF while transposing\n",
    "        tempDF = pd.DataFrame(tempS.to_dict(), index=[0])\n",
    "        ##Clean DF\n",
    "        tempDF['N'] = tempDF['N'].astype('int64')\n",
    "        tempDF['DoF'] = tempDF['DoF'].astype('int64')\n",
    "        tempL = tempDF.loc[:, tempDF.columns.str.contains('_nM.*ed')].columns.tolist()\n",
    "        for col_n in tempL:\n",
    "            tempDF[col_n] = tempDF[col_n].astype('int64')\n",
    "        tempDF['BMIclass'] = bmi_class\n",
    "        \n",
    "        tempD2[bmi_class] = tempDF\n",
    "    \n",
    "    #Clean the results (pd.DataFrame) across BMI classes\n",
    "    tempDF = pd.concat(list(tempD2.values()), axis=0)\n",
    "    tempDF['Metric'] = metric.replace('_tf', '')\n",
    "    \n",
    "    #P-value adjustment (across BMI classes and bBMIs within the metric) by using Benjamini–Hochberg method\n",
    "    tempL = tempDF.loc[:, tempDF.columns.str.contains('_Pval')].columns.tolist()\n",
    "    tempDF4 = tempDF.reset_index().melt(var_name='bBMI', value_name='Pval', value_vars=tempL,\n",
    "                                        id_vars=['BMIclass', 'Metric'])\n",
    "    tempDF4['AdjPval'] = multi.multipletests(tempDF4['Pval'], alpha=0.05, method='fdr_bh',\n",
    "                                             is_sorted=False, returnsorted=False)[1]\n",
    "    tempDF4 = tempDF4.pivot(index=['BMIclass', 'Metric'], columns='bBMI', values='AdjPval')\n",
    "    tempDF4.columns = tempDF4.columns.str.replace('_Pval', '_AdjPval_within1')\n",
    "    ##Replace the dummy values with the adjusted p-values\n",
    "    tempL = [(bmi_class, metric.replace('_tf', '')) for bmi_class in tempL1]\n",
    "    tempDF4 = tempDF4.loc[tempL]#Sort just in case\n",
    "    tempL = tempDF.loc[:, tempDF.columns.str.contains('_AdjPval_within1')].columns.tolist()\n",
    "    for col_n in tempL:\n",
    "        tempDF[col_n] = tempDF4[col_n].tolist()\n",
    "    \n",
    "    tempD1[metric] = tempDF\n",
    "t_elapsed = time.time() - t_start\n",
    "print('Elapsed time for',\n",
    "      len(tempDF1.columns)*len(tempL1)*len(tempL2), 'OLS linear regressions (',\n",
    "      len(tempDF1.columns), 'metrics x',\n",
    "      len(tempL1), 'BMI classes x',\n",
    "      len(tempL2), 'bBMIs):',\n",
    "      round(t_elapsed//60), 'min', round(t_elapsed%60, 1), 'sec')\n",
    "\n",
    "#Clean the results (pd.DataFrame) across metrics\n",
    "tempDF = pd.concat(list(tempD1.values()), axis=0)\n",
    "##Clean the column order by setting index\n",
    "tempDF = tempDF.set_index(['BMIclass', 'Metric'])\n",
    "\n",
    "#P-value adjustment (across BMI classes and metrics within each bBMI) by using Benjamini–Hochberg method\n",
    "for bbmi in tempL2:\n",
    "    tempDF['vs'+bbmi+'class_AdjPval_within2'] = multi.multipletests(tempDF['vs'+bbmi+'class_Pval'],\n",
    "                                                                    alpha=0.05, method='fdr_bh',\n",
    "                                                                    is_sorted=False, returnsorted=False)[1]\n",
    "\n",
    "#P-value adjustment (across all tests) by using Benjamini–Hochberg method\n",
    "tempL = tempDF.loc[:, tempDF.columns.str.contains('_Pval')].columns.tolist()\n",
    "tempDF4 = tempDF.reset_index().melt(var_name='bBMI', value_name='Pval', value_vars=tempL,\n",
    "                                    id_vars=['BMIclass', 'Metric'])\n",
    "tempDF4['AdjPval'] = multi.multipletests(tempDF4['Pval'], alpha=0.05, method='fdr_bh',\n",
    "                                         is_sorted=False, returnsorted=False)[1]\n",
    "tempDF4 = tempDF4.pivot(index=['BMIclass', 'Metric'], columns='bBMI', values='AdjPval')\n",
    "tempDF4.columns = tempDF4.columns.str.replace('_Pval', '_AdjPval_all')\n",
    "##Replace the dummy values with the adjusted p-values\n",
    "tempL = tempDF.loc[:, tempDF.columns.str.contains('_AdjPval_all')].columns.tolist()\n",
    "for col_n in tempL:\n",
    "    tempDF[col_n] = tempDF4[col_n]\n",
    "\n",
    "tempDF = tempDF.sort_index(axis=0, ascending=True, key=lambda x:x.str.lower())\n",
    "display(tempDF)\n",
    "\n",
    "#Save\n",
    "fileDir = './ExportData/'\n",
    "ipynbName = '220905_Multiomics-BMI-NatMed1stRevision_Microbiome-DeltaBMI-wenceslaus_'\n",
    "fileName = 'result-summary.tsv'\n",
    "tempDF.to_csv(fileDir+ipynbName+fileName, sep='\\t', index=True)\n",
    "\n",
    "resDF = tempDF"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2-2. Significantly different metrics by misclassification"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempDF = resDF\n",
    "tempL1 = ['Normal', 'Obese']\n",
    "tempL2 = ['MetBMI', 'ProtBMI', 'ChemBMI', 'CombiBMI']\n",
    "\n",
    "#Significantly different metrics\n",
    "print('Significantly different metrics by misclassification (FDR < 0.05)')\n",
    "tempD = {}\n",
    "for bmi_class in tempL1:\n",
    "    tempDF1 = tempDF.loc[bmi_class]#MultiIndex\n",
    "    print(' - '+bmi_class)\n",
    "    tempL = []\n",
    "    for bbmi in tempL2:\n",
    "        #Extact significant metrics\n",
    "        tempDF2 = tempDF1.loc[tempDF1['vs'+bbmi+'class_AdjPval_all']<0.05]\n",
    "        tempL.append(tempDF2.index.tolist())\n",
    "        \n",
    "        #Cf. Adjustment across BMI classes and bBMIs within the metric\n",
    "        tempDF3 = tempDF1.loc[tempDF1['vs'+bbmi+'class_AdjPval_within1']<0.05]\n",
    "        #Cf. Adjustment across BMI classes and metrics within each bBMI\n",
    "        tempDF4 = tempDF1.loc[tempDF1['vs'+bbmi+'class_AdjPval_within2']<0.05]\n",
    "        #Cf. Nominal P-value\n",
    "        tempDF5 = tempDF1.loc[tempDF1['vs'+bbmi+'class_Pval']<0.05]\n",
    "        #tempL.append(tempDF5.index.tolist())\n",
    "        \n",
    "        print('   - '+bbmi+':', len(tempDF2),\n",
    "              '(cf. within-adjustment 1:', len(tempDF3),\n",
    "              ', within-adjustment 2:', len(tempDF4),\n",
    "              ', nominal P-value:', len(tempDF5), ')')\n",
    "    #Flatten and drop multiplicates\n",
    "    tempL = list(set(row_n for sublist in tempL for row_n in sublist))\n",
    "    tempDF1 = tempDF1.loc[tempL]\n",
    "    \n",
    "    print('   - Union:', len(tempL))\n",
    "    tempD[bmi_class] = tempL\n",
    "    \n",
    "    #Clean (just for the display in Jupyter notebook)\n",
    "    tempL3 = tempDF1.loc[:, tempDF1.columns.str.contains('Bcoef$')]\n",
    "    tempL4 = tempDF1.loc[:, tempDF1.columns.str.contains('Pval')]\n",
    "    tempL = [col_n for sublist in [tempL3, tempL4] for col_n in sublist]\n",
    "    tempDF1 = tempDF1[tempL]\n",
    "    tempDF1 = tempDF1.sort_values(by='vsMetBMIclass_Pval', ascending=True)\n",
    "    display(tempDF1)\n",
    "\n",
    "#Flatten and drop multiplicates\n",
    "tempL = list(set(row_n for sublist in tempD.values() for row_n in sublist))\n",
    "print(' - Union:', len(tempL))\n",
    "tempL = [(bmi_class, metric) for bmi_class in tempL1 for metric in tempL]\n",
    "tempDF1 = tempDF.loc[tempL]\n",
    "#Clean (just for the display in Jupyter notebook)\n",
    "tempL3 = tempDF1.loc[:, tempDF1.columns.str.contains('Bcoef$')]\n",
    "tempL4 = tempDF1.loc[:, tempDF1.columns.str.contains('Pval')]\n",
    "tempL = [col_n for sublist in [tempL3, tempL4] for col_n in sublist]\n",
    "tempDF1 = tempDF1[tempL]\n",
    "tempDF1 = tempDF1.sort_values(by='vsMetBMIclass_Pval', ascending=True)\n",
    "display(tempDF1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2-3. Visualization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempD1 = {'Observed_tf':'Observed ASVs',\n",
    "          'Shannon_tf':'Shannon\\'s index',\n",
    "          'Chao1_tf':'Chao1 diversity'}\n",
    "tempD2 = {'Underweight':'blue', 'Normal':'green', 'Overweight':'orange', 'Obese':'red'}\n",
    "tempD3 = {'MetBMI':'b', 'ProtBMI':'r', 'ChemBMI':'g', 'CombiBMI':'m'}\n",
    "tempD4 = {'Observed_tf':'tab:blue',\n",
    "          'Shannon_tf':'tab:blue',\n",
    "          'Chao1_tf':'tab:blue'}\n",
    "tempL1 = ['Normal', 'Obese']\n",
    "tempDF1 = divDF\n",
    "tempDF2 = bmiDF\n",
    "tempDF3 = resDF\n",
    "\n",
    "for metric_i, metric in enumerate(tempD1.keys()):\n",
    "    #Prepare DF\n",
    "    tempS = tempDF1[metric]\n",
    "    tempDF = pd.merge(tempS, tempDF2, left_index=True, right_index=True, how='left')\n",
    "    ##Drop NaN in the metric values\n",
    "    #tempDF = tempDF.dropna()\n",
    "    ##Z-score transformation (based on the whole distribution)\n",
    "    tempDF4 = tempDF.select_dtypes(include=[np.number])\n",
    "    scaler = StandardScaler(copy=True, with_mean=True, with_std=True)\n",
    "    tempA = scaler.fit_transform(tempDF4)#Column direction\n",
    "    tempDF4 = pd.DataFrame(data=tempA, index=tempDF4.index, columns=tempDF4.columns)\n",
    "    ###Recover the categorical variables\n",
    "    tempDF5 = tempDF.select_dtypes(exclude=[np.number])\n",
    "    tempDF = pd.merge(tempDF4, tempDF5, left_index=True, right_index=True, how='left')\n",
    "    \n",
    "    #Check sample size\n",
    "    print(tempD1[metric])\n",
    "    print('N (total):', len(tempDF))\n",
    "    print(' - BMI class:', tempDF['BMI_class'].value_counts().sort_index(ascending=True).to_dict())\n",
    "    for bmi_class in tempL1:\n",
    "        print('   - '+bmi_class+' BMI class')\n",
    "        for bbmi in tempD3.keys():\n",
    "            tempS = tempDF['vs_'+bbmi+'_class'].loc[tempDF['BMI_class']==bmi_class]\n",
    "            print('     - vs. '+bbmi+' class:',\n",
    "                  tempS.value_counts().sort_index(ascending=True).to_dict())\n",
    "    \n",
    "    #Visualization\n",
    "    sns.set(style='ticks', font='Arial', context='talk')\n",
    "    fig, axes = plt.subplots(nrows=1, ncols=1+len(tempD3),\n",
    "                             figsize=(9, 3), sharex=False, sharey=True,\n",
    "                             gridspec_kw={'width_ratios':[1, 1, 1, 1, 1]})\n",
    "    axis_ymin = -3.5\n",
    "    axis_ymax = 3.75\n",
    "    ymin = -3\n",
    "    ymax = 3\n",
    "    yinter = 1\n",
    "    margin = 0.49\n",
    "    #Set shared axis range\n",
    "    plt.setp(axes, ylim=(axis_ymin, axis_ymax), yticks=np.arange(ymin, ymax+yinter/10, yinter))\n",
    "    tempL = []#For legend handles and labels\n",
    "    for ax_i, ax in enumerate(axes.flat):\n",
    "        if ax_i==0:#Overall\n",
    "            sns.boxplot(data=tempDF, y=metric, x='BMI_class', order=tempD2.keys(),\n",
    "                        hue='BMI_class', hue_order=tempD2.keys(), dodge=False, palette=tempD2,\n",
    "                        showfliers=False,#flierprops={'marker':'o', 'markerfacecolor':'gray', 'alpha':0.4},\n",
    "                        showcaps=True, notch=True, ax=ax)\n",
    "            tempL.append(ax.get_legend_handles_labels())\n",
    "        else:#Misclassification\n",
    "            bbmi = list(tempD3.keys())[ax_i-1]\n",
    "            tempDF4 = tempDF.loc[tempDF['BMI_class'].isin(tempL1)]\n",
    "            tempD = {'Matched':'0.8', 'Mismatched':tempD3[bbmi]}\n",
    "            sns.boxplot(data=tempDF4, y=metric, x='BMI_class', order=tempL1,\n",
    "                        hue='vs_'+bbmi+'_class', hue_order=tempD.keys(), dodge=True, palette=tempD,\n",
    "                        showfliers=False,#flierprops={'marker':'o', 'markerfacecolor':'gray', 'alpha':0.4},\n",
    "                        showcaps=True, notch=True, ax=ax)\n",
    "            tempL.append(ax.get_legend_handles_labels())\n",
    "        #Axis settings\n",
    "        if ax_i==0:\n",
    "            plt.setp(ax, xlim=(0-margin, len(tempD2)-1+margin))#To eliminate excess white space\n",
    "            plt.setp(ax, xlabel='', ylabel='Tarnsformed value [a.u.]\\n('+r'$Z$'+'-score)')\n",
    "        else:\n",
    "            plt.setp(ax, xlim=(0-margin, len(tempD)-1+margin))#To eliminate excess white space\n",
    "            plt.setp(ax.get_yticklabels(), visible=False)\n",
    "            plt.setp(ax, xlabel='', ylabel='')\n",
    "        #ax.grid(axis='y', linestyle='--', color='black')\n",
    "        sns.despine()\n",
    "        plt.setp(ax.get_xticklabels(), rotation=70,\n",
    "                 horizontalalignment='right', verticalalignment='center', rotation_mode='anchor')\n",
    "        #Annotation\n",
    "        lines = ax.get_lines()#Line2D: [[Q1, Q1-1.5IQR], [Q3, Q3+1.5IQR], [Q1, Q1], [Q3, Q3], [Med, Med], [flier]]\n",
    "        if ax_i!=0:\n",
    "            #P-value annotation\n",
    "            lines_unit = 5 + int(False)#showfliers=False\n",
    "            for class_i in range(len(tempL1)):\n",
    "                #Matched\n",
    "                whisker_0 = lines[class_i*lines_unit*len(tempD) + lines_unit*0 + 1]\n",
    "                xcoord_0 = whisker_0._x[1]#Q3+1.5IQR\n",
    "                ycoord_0 = whisker_0._y[1]#Q3+1.5IQR\n",
    "                #Mismatched\n",
    "                whisker_1 = lines[class_i*lines_unit*len(tempD) + lines_unit*1 + 1]\n",
    "                xcoord_1 = whisker_1._x[1]#Q3+1.5IQR\n",
    "                ycoord_1 = whisker_1._y[1]#Q3+1.5IQR\n",
    "                #Standard point for annotation\n",
    "                xcoord = (xcoord_0+xcoord_1)/2\n",
    "                ycoord = max(ycoord_0, ycoord_1)\n",
    "                #Add annotation lines\n",
    "                aline_offset = yinter/5\n",
    "                aline_length = yinter/5 + aline_offset/2\n",
    "                ax.plot([xcoord_0, xcoord_0, xcoord_1, xcoord_1],\n",
    "                        [ycoord+aline_offset, ycoord+aline_length, ycoord+aline_length, ycoord+aline_offset],\n",
    "                        lw=1.5, c='k')\n",
    "                #Retrieve P-value\n",
    "                bmi_class = tempL1[class_i]\n",
    "                pval = tempDF3.loc[(bmi_class, metric.replace('_tf', '')), 'vs'+bbmi+'class_AdjPval_all']\n",
    "                if pval<0.001:\n",
    "                    label = '***'\n",
    "                elif pval<0.01:\n",
    "                    label = '**'\n",
    "                elif pval<0.05:\n",
    "                    label = '*'\n",
    "                else:\n",
    "                    pval_text = str(Decimal(pval).quantize(Decimal('0.01'), rounding=ROUND_HALF_UP))\n",
    "                    label = r'$P$'+' = '+pval_text\n",
    "                #Add annotation text\n",
    "                if label in ['***', '**', '*']:\n",
    "                    text_offset = yinter/12\n",
    "                    text_size = 'medium'\n",
    "                else:\n",
    "                    text_offset = yinter/3\n",
    "                    text_size = 'x-small'\n",
    "                ax.annotate(label, xy=(xcoord, ycoord+text_offset),\n",
    "                            horizontalalignment='center', verticalalignment='bottom',\n",
    "                            fontsize=text_size, color='k')\n",
    "        #Facet settings\n",
    "        if ax_i==0:\n",
    "            ax.set_title('Overall', {'fontsize':'medium'})\n",
    "        else:\n",
    "            ax.set_title(bbmi, {'fontsize':'medium'})\n",
    "            xoff = 0.025\n",
    "            yoff = 0.01\n",
    "            rect = plt.Rectangle((xoff, 1+yoff), 1-xoff, 0.15,#Manual adjustment\n",
    "                                 transform=ax.transAxes, facecolor=tempD3[bbmi], alpha=0.3,\n",
    "                                 clip_on=False, linewidth=0, zorder=0.5)\n",
    "            ax.add_patch(rect)\n",
    "        #Change the default boxplot settings\n",
    "        for line in lines:\n",
    "            line.set_color('k')\n",
    "        for box in ax.artists:\n",
    "            box.set_edgecolor('k')\n",
    "        #Remove the default legend\n",
    "        ax.get_legend().remove()\n",
    "        #Save ax position for figure title\n",
    "        if ax_i==0:\n",
    "            ax_pos_l = ax.get_position().bounds\n",
    "        elif ax_i==len(tempD3):\n",
    "            ax_pos_r = ax.get_position().bounds\n",
    "    #Add legend manually\n",
    "    for legend_i in range(len(tempL)):\n",
    "        h, l = tempL[legend_i]\n",
    "        if legend_i==0:\n",
    "            title_text = 'BMI class (overall)'\n",
    "            position = (0.075, -0.3)\n",
    "        else:\n",
    "            bbmi =list(tempD3.keys())[legend_i-1]\n",
    "            title_text = 'vs. '+bbmi+' class'\n",
    "            position = (0.35+0.25*((legend_i-1)//2), -0.3-0.35*((legend_i-1)%2))\n",
    "        legend = fig.legend(handles=h, labels=l, fontsize='medium',\n",
    "                            title=title_text, title_fontsize='medium',\n",
    "                            bbox_to_anchor=position, loc='upper left',\n",
    "                            labelspacing=0.25, handletextpad=0.5,\n",
    "                            borderaxespad=1, frameon=False)\n",
    "        plt.gca().add_artist(legend)\n",
    "    #Add figure title\n",
    "    xcoord_0 = ax_pos_l[0]\n",
    "    xcoord_1 = ax_pos_r[0]+ax_pos_r[2]\n",
    "    ycoord = ax_pos_l[1]+ax_pos_l[3]\n",
    "    yoff = 0.15\n",
    "    if tempD4[metric] in ['tab:red', 'tab:blue']:\n",
    "        text_color = 'white'\n",
    "    else:\n",
    "        text_color = 'black'\n",
    "    fig.suptitle(tempD1[metric], x=(xcoord_0+xcoord_1)/2, y=ycoord+yoff,\n",
    "                 fontsize='large', fontweight='bold', color=text_color,\n",
    "                 horizontalalignment='center', verticalalignment='bottom')\n",
    "    yoff = yoff - 0.015#Manual adjustment\n",
    "    rect = plt.Rectangle((xcoord_0, ycoord+yoff), xcoord_1-xcoord_0, 0.125,#Manual adjustment\n",
    "                         transform=fig.transFigure, facecolor=tempD4[metric],\n",
    "                         clip_on=False, linewidth=0, zorder=0)\n",
    "    fig.patches.extend([rect])\n",
    "    ##Save\n",
    "    fileDir = './ExportFigures/'\n",
    "    ipynbName = '220905_Multiomics-BMI-NatMed1stRevision_Microbiome-DeltaBMI-wenceslaus_'\n",
    "    fileName = metric.replace('_tf', '')+'.tif'\n",
    "    plt.gcf().savefig(fileDir+ipynbName+fileName, dpi=300, bbox_inches='tight', pad_inches=0.04,\n",
    "                      pil_kwargs={'compression':'tiff_lzw'})\n",
    "    plt.show()\n",
    "    print('')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# — End of this notebook —"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "arivale-py3 - Python",
   "language": "python",
   "name": "conda-env-arivale-py3-py"
  },
  "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.9.7"
  },
  "toc-showcode": false,
  "toc-showmarkdowntxt": false,
  "toc-showtags": false,
  "widgets": {
   "application/vnd.jupyter.widget-state+json": {
    "state": {},
    "version_major": 2,
    "version_minor": 0
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}