Code08_DeltaBMI-misclassification_clustering.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Multiomics BMI Paper — ∆BMI-based Misclassification and Hierarchical Clustering"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***by Kengo Watanabe***  \n",
    "\n",
    "This Jupyter Notebook (with Python 3 kernel) visualized ∆BMI correlations across omics categories and ∆BMI-based misclassification (in the baseline Arivale cohort). Also, this notebook performed hierarchical clustering to reveal the underlying differences in ∆BMI misclassification, which was removed from the final paper but included in the original preprint.  \n",
    "\n",
    "Input files:  \n",
    "* Arivale baseline BMI predictions: 220801_Multiomics-BMI-NatMed1stRevision_BMI-baseline-LASSO_\\[MetBMI/ProtBMI/ChemBMI/CombiBMI\\]-BothSex.tsv  \n",
    "* Arivale baseline covariates and blood omics (preprocessed): 210104_Biological-BMI-paper_RF-imputation_baseline-\\[metDF/protDF/combiDF\\]-with-RF-imputation.tsv  \n",
    "* OLS linear regressions of the retained variables on BMI: 220802_Multiomics-BMI-NatMed1stRevision_BMI-LASSO-bcoef_\\[MetBMI/ProtBMI\\]-BothSex-OLS.tsv  \n",
    "\n",
    "Output figures and tables:  \n",
    "* Figure 3a, 3c  \n",
    "* Tables for Supplementary Data 10  \n",
    "* Intermediate table for other notebooks (biological BMI summary)  \n",
    "* (Figure 3d, 3e and Supplementary Figure 5b, 5c in the original preprint)  \n",
    "\n",
    "Original notebook (memo for my future tracing):  \n",
    "* dalek:\\[JupyterLab HOME\\]/220621_Multiomics-BMI-NatMedRevision/220803_Multiomics-BMI-NatMed1stRevision_DeltaBMI-misclassification.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            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        py39hf3d152e_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",
      "matplotlib-venn           0.11.7             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.1o               h166bdaf_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",
      "venn                      0.1.3                    pypi_0    pypi\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 decimal import Decimal, ROUND_HALF_UP\n",
    "import sys\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from scipy.cluster import hierarchy\n",
    "from statsmodels.stats import multitest as multi\n",
    "import matplotlib.patches as mpatches\n",
    "\n",
    "!conda list"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Obesity classification"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1-1. Clean the biological BMI dataframes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Import and merge BMI and biological BMI\n",
    "fileDir = './ExportData/'\n",
    "ipynbName = '220801_Multiomics-BMI-NatMed1stRevision_BMI-baseline-LASSO_'\n",
    "tempL1 = ['log_BaseBMI', 'BaseBMI', 'Testing']\n",
    "tempL2 = ['MetBMI', 'ProtBMI', 'ChemBMI', 'CombiBMI']\n",
    "for bbmi in tempL2:\n",
    "    fileName = bbmi+'-BothSex.tsv'\n",
    "    tempDF1 = pd.read_csv(fileDir+ipynbName+fileName, sep='\\t', dtype={'public_client_id':str})\n",
    "    tempDF1 = tempDF1.set_index('public_client_id')\n",
    "    #Put aside common part\n",
    "    if bbmi==tempL2[0]:\n",
    "        tempDF = tempDF1[tempL1]\n",
    "    tempDF1 = tempDF1.drop(columns=tempL1)\n",
    "    tempDF = pd.merge(tempDF, tempDF1, left_index=True, right_index=True, how='inner')\n",
    "\n",
    "display(tempDF)\n",
    "tempDF1 = tempDF.describe(include=[np.number])\n",
    "tempDF1.loc['Skewness'] = stats.skew(tempDF.select_dtypes(include=[np.number]))\n",
    "display(tempDF1)\n",
    "\n",
    "tempD = {'BMI':'k', 'MetBMI':'b', 'ProtBMI':'r', 'ChemBMI':'g', 'CombiBMI':'m'}\n",
    "for scale in ['log', 'raw']:\n",
    "    sns.set(style='ticks', font='Arial', context='talk')\n",
    "    plt.figure(figsize=(4, 3))\n",
    "    for bbmi in tempD.keys():\n",
    "        if scale=='raw':\n",
    "            sns.distplot(tempDF['Base'+bbmi], label=bbmi, color=tempD[bbmi])\n",
    "        elif scale=='log':\n",
    "            sns.distplot(tempDF['log_Base'+bbmi], label=bbmi, color=tempD[bbmi])\n",
    "    sns.despine()\n",
    "    plt.ylabel('Density')\n",
    "    if scale=='raw':\n",
    "        plt.xlabel('BMI [kg/m'+r'$^2$'+']')\n",
    "    elif scale=='log':\n",
    "        plt.xlabel('BMI [kg/m'+r'$^2$'+'] (log-scale)')\n",
    "    plt.legend(bbox_to_anchor=(1, 0.5), loc='center left', borderaxespad=1)\n",
    "    plt.show()\n",
    "\n",
    "bmiDF = tempDF"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1-2. Obesity classification"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Obesity classification\n",
    "tempD = {'BMI':'k', 'MetBMI':'b', 'ProtBMI':'r', 'ChemBMI':'g', 'CombiBMI':'m'}\n",
    "for bmi in tempD.keys():\n",
    "    tempL = []\n",
    "    for value in bmiDF['Base'+bmi].tolist():\n",
    "        if np.isnan(value):\n",
    "            tempL.append('NotCalculated')\n",
    "        elif value < 18.5:\n",
    "            tempL.append('Underweight')\n",
    "        elif value < 25:\n",
    "            tempL.append('Normal')\n",
    "        elif value < 30:\n",
    "            tempL.append('Overweight')\n",
    "        elif value >= 30:\n",
    "            tempL.append('Obese')\n",
    "        else:#Just in case\n",
    "            tempL.append('Error?')\n",
    "    bmiDF['Base'+bmi+'_class'] = tempL\n",
    "    #Confirmation\n",
    "    print('Base'+bmi+'_class:')\n",
    "    tempS1 = bmiDF['Base'+bmi+'_class'].value_counts()\n",
    "    tempDF = pd.DataFrame({'Count':tempS1, 'Percentage':tempS1/len(bmiDF)*100})\n",
    "    display(tempDF)\n",
    "    print('')\n",
    "\n",
    "#Add the covariates info to make summary dataframe\n",
    "fileDir = '../210104_Biological-BMI-paper/ExportData/'\n",
    "ipynbName = '210104_Biological-BMI-paper_RF-imputation_'\n",
    "fileName = 'baseline-combiDF-with-RF-imputation.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",
    "tempL = ['Sex', 'BaseAge', 'PC1', 'PC2', 'PC3', 'PC4', 'PC5']#Without Race in this study\n",
    "bmiDF = pd.merge(bmiDF, tempDF[tempL], left_index=True, right_index=True, how='inner')\n",
    "display(bmiDF)\n",
    "\n",
    "#Save\n",
    "fileDir = './ExportData/'\n",
    "ipynbName = '220803_Multiomics-BMI-NatMed1stRevision_DeltaBMI-misclassification_'\n",
    "fileName = 'biologicalBMI-baseline-summary-BothSex.tsv'\n",
    "bmiDF.to_csv(fileDir+ipynbName+fileName, sep='\\t', index=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1-3. Misclassification"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Calculate misclassification rate based on each biological BMI class\n",
    "tempL = ['MetBMI', 'ProtBMI', 'ChemBMI', 'CombiBMI']\n",
    "misclassDF = pd.DataFrame(index=['Overall_count', 'Overall [%]',\n",
    "                                 'Underweight_count', 'Underweight [%]',\n",
    "                                 'Normal_count', 'Normal [%]',\n",
    "                                 'Overweight_count', 'Overweight [%]',\n",
    "                                 'Obese_count', 'Obese [%]'])\n",
    "for bbmi in tempL:\n",
    "    tempDF = bmiDF.loc[bmiDF['Base'+bbmi+'_class']!='NotCalculated']#just in case\n",
    "    counter0 = 0\n",
    "    counter1 = 0\n",
    "    counter2 = 0\n",
    "    counter3 = 0\n",
    "    counter4 = 0\n",
    "    for row_n in tempDF.index.tolist():\n",
    "        bmi_class = tempDF.loc[row_n, 'BaseBMI_class']\n",
    "        bbmi_class = tempDF.loc[row_n, 'Base'+bbmi+'_class']\n",
    "        if bmi_class!=bbmi_class:\n",
    "            counter0 += 1#Overall\n",
    "            if bmi_class=='Underweight':\n",
    "                counter1 += 1\n",
    "            elif bmi_class=='Normal':\n",
    "                counter2 += 1\n",
    "            elif bmi_class=='Overweight':\n",
    "                counter3 += 1\n",
    "            elif bmi_class=='Obese':\n",
    "                counter4 += 1\n",
    "            else:#Just in case\n",
    "                print('Error?')\n",
    "    tempS = pd.Series([counter0, counter0/len(tempDF)*100,\n",
    "                       counter1, counter1/len(tempDF.loc[tempDF['BaseBMI_class']=='Underweight'])*100,\n",
    "                       counter2, counter2/len(tempDF.loc[tempDF['BaseBMI_class']=='Normal'])*100,\n",
    "                       counter3, counter3/len(tempDF.loc[tempDF['BaseBMI_class']=='Overweight'])*100,\n",
    "                       counter4, counter4/len(tempDF.loc[tempDF['BaseBMI_class']=='Obese'])*100],\n",
    "                      index=['Overall_count', 'Overall [%]',\n",
    "                             'Underweight_count', 'Underweight [%]',\n",
    "                             'Normal_count', 'Normal [%]',\n",
    "                             'Overweight_count', 'Overweight [%]',\n",
    "                             'Obese_count', 'Obese [%]'],\n",
    "                      name='vs. '+bbmi)\n",
    "    misclassDF = pd.concat([misclassDF, tempS], axis=1)\n",
    "display(misclassDF)\n",
    "\n",
    "#Plot (without underweight)\n",
    "tempD1 = {'vs. MetBMI':'vs. MetBMI class', 'vs. ProtBMI':'vs. ProtBMI class',\n",
    "          'vs. ChemBMI':'vs. ChemBMI class', 'vs. CombiBMI':'vs. CombiBMI class'}\n",
    "tempD2 = {'vs. MetBMI class':'b', 'vs. ProtBMI class':'r',\n",
    "          'vs. ChemBMI class':'g', 'vs. CombiBMI class':'m'}\n",
    "tempDF = misclassDF.loc[['Overall [%]', 'Normal [%]', 'Overweight [%]', 'Obese [%]']]\n",
    "tempDF.index = tempDF.index.str.replace(' \\[%\\]', '')\n",
    "tempDF.columns = tempDF.columns.map(tempD1)\n",
    "tempDF = tempDF.reset_index().melt(var_name='Category', value_name='Misclassification', id_vars='index')\n",
    "sns.set(style='ticks', font='Arial', context='talk')\n",
    "plt.figure(figsize=(5.5, 3))\n",
    "sns.barplot(data=tempDF, x='index', y='Misclassification',\n",
    "            hue='Category', dodge=True, palette=tempD2, edgecolor='black')\n",
    "sns.despine()\n",
    "plt.yticks(np.arange(0, 51, 10))\n",
    "plt.axvline(x=(0+1)/2, **{'linestyle':'--', 'color':'k'})\n",
    "##Add reference range\n",
    "plt.axhspan(ymin=28, ymax=48, facecolor='orange', alpha=0.2, zorder=0)\n",
    "plt.axhline(y=28, **{'linestyle':'-', 'color':'orange', 'zorder':0})\n",
    "plt.axhline(y=48, **{'linestyle':'-', 'color':'orange', 'zorder':0})\n",
    "plt.ylabel('Misclassification [%]')\n",
    "plt.xlabel('BMI class')\n",
    "#plt.legend(bbox_to_anchor=(0.5, 1), loc='lower center', borderaxespad=0.5,\n",
    "#           ncol=2, columnspacing=1, handletextpad=0.5)\n",
    "plt.legend(bbox_to_anchor=(1, 0.5), loc='center left', borderaxespad=0.5)\n",
    "##Save\n",
    "fileDir = './ExportFigures/'\n",
    "ipynbName = '220803_Multiomics-BMI-NatMed1stRevision_DeltaBMI-misclassification_'\n",
    "fileName = 'misclassification-rate.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()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1-4. Calculate ∆BMI"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> In contrast to biological age, the difference between measured and predicted values are dependent on the absoculte BMI value. Therefore, ∆BMI is defined with the rate of difference: (bBMI - BMI) / BMI.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Calculate the rate of difference\n",
    "tempL = ['MetBMI', 'ProtBMI', 'ChemBMI', 'CombiBMI']\n",
    "for bbmi in tempL:\n",
    "    bmiDF[bbmi+'–BMI'] = (bmiDF['Base'+bbmi] - bmiDF['BaseBMI']) / bmiDF['BaseBMI'] * 100\n",
    "\n",
    "tempD = {'MetBMI–BMI':'b', 'ProtBMI–BMI':'r', 'ChemBMI–BMI':'g', 'CombiBMI–BMI':'m'}\n",
    "tempDF = bmiDF[list(tempD.keys())]\n",
    "tempDF1 = tempDF.describe()\n",
    "tempDF1.loc['Skewness'] = stats.skew(tempDF)\n",
    "display(tempDF1)\n",
    "\n",
    "sns.set(style='ticks', font='Arial', context='talk')\n",
    "plt.figure(figsize=(4, 3))\n",
    "for col_n in tempD.keys():\n",
    "    sns.distplot(tempDF[col_n], label=col_n, color=tempD[col_n])\n",
    "sns.despine()\n",
    "plt.ylabel('Density')\n",
    "plt.xlabel(r'$\\Delta$'+'BMI [% BMI]')\n",
    "plt.legend(bbox_to_anchor=(1, 0.5), loc='center left', borderaxespad=1)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1-5. Difference in ∆BMI between omics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> In this version, P-values from Pearson's correlation tests are adjusted across all comparisons.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Prepare DF\n",
    "tempL = ['MetBMI–BMI', 'ProtBMI–BMI', 'ChemBMI–BMI', 'CombiBMI–BMI', 'BaseBMI_class']\n",
    "tempDF = bmiDF[tempL]\n",
    "tempDF.columns = tempDF.columns.str.replace('–BMI', '')\n",
    "\n",
    "#Check correlation matrix and extract lower triangle matrix\n",
    "tempDF1 = tempDF.select_dtypes(include=[np.number])\n",
    "tempDF1 = tempDF1.corr(method='pearson')\n",
    "print('Pearson\\'s r:')\n",
    "display(tempDF1)\n",
    "tempDF1 = tempDF1.where(np.tril(np.ones(tempDF1.shape), k=-1).astype(np.bool), other=np.nan)\n",
    "tempDF1.index.rename('Variable1', inplace=True)\n",
    "tempDF1 = tempDF1.reset_index().melt(var_name='Variable2', value_name='Pearson_r', id_vars=['Variable1'])\n",
    "tempDF1 = tempDF1.dropna()\n",
    "\n",
    "#Statistical tests\n",
    "tempDF2 = pd.DataFrame(columns=['Xvar', 'Yvar', 'N', 'DoF', 'Pearson_r', 'Pval'])\n",
    "for row_i in range(len(tempDF1)):\n",
    "    xvar = tempDF1['Variable2'].iloc[row_i]\n",
    "    yvar = tempDF1['Variable1'].iloc[row_i]\n",
    "    #Pearson's correlation\n",
    "    pearson_r, pval = stats.pearsonr(tempDF[xvar], tempDF[yvar])\n",
    "    size = len(tempDF)\n",
    "    dof = size - 2\n",
    "    tempDF2.loc[xvar+'-vs-'+yvar] = [xvar, yvar, size, dof, pearson_r, pval]\n",
    "##P-value adjustment by using Benjamini–Hochberg method\n",
    "tempDF2['AdjPval'] = multi.multipletests(tempDF2['Pval'], alpha=0.05, method='fdr_bh',\n",
    "                                         is_sorted=False, returnsorted=False)[1]\n",
    "tempDF2.index.rename('ComparisonLabel', inplace=True)\n",
    "display(tempDF2)\n",
    "##Save\n",
    "fileDir = './ExportData/'\n",
    "ipynbName = '220803_Multiomics-BMI-NatMed1stRevision_DeltaBMI-misclassification_'\n",
    "fileName = 'bBMI-BMI-difference.tsv'\n",
    "tempDF2.to_csv(fileDir+ipynbName+fileName, sep='\\t', index=True)\n",
    "\n",
    "#Visualization\n",
    "tempD = {'Underweight':'blue', 'Normal':'green', 'Overweight':'orange', 'Obese':'red'}\n",
    "sns.set(style='ticks', font='Arial', context='talk')\n",
    "p = sns.PairGrid(tempDF, hue='BaseBMI_class', hue_order=list(tempD.keys()), palette=tempD,\n",
    "                 height=2, aspect=1, layout_pad=0.0)\n",
    "p.map_lower(sns.scatterplot, edgecolor='0.3', alpha=0.5, s=25)\n",
    "p.map_diag(sns.distplot, axlabel=False, kde_kws={'alpha':0.8}, hist_kws={'edgecolor':'white', 'alpha':0.5})\n",
    "for i, j in zip(*np.triu_indices_from(p.axes, 1)):\n",
    "    p.axes[i, j].set_visible(False)\n",
    "for i, j in zip(*np.tril_indices_from(p.axes, 0)):\n",
    "    p.axes[i, j].set(xlim=(-40, 85), xticks=np.arange(-25, 75.1, 25),\n",
    "                     ylim=(-40, 85), yticks=np.arange(-25, 75.1, 25))\n",
    "for i, j in zip(*np.tril_indices_from(p.axes, -1)):\n",
    "    p.axes[i, j].grid(axis='both', linestyle='--', color='gray', alpha=0.3)\n",
    "    #Annotate Pearson's correlation\n",
    "    xvar = tempDF.columns.tolist()[j]\n",
    "    yvar = tempDF.columns.tolist()[i]\n",
    "    pearson_r = tempDF2['Pearson_r'].loc[xvar+'-vs-'+yvar]\n",
    "    r_text = str(Decimal(str(pearson_r)).quantize(Decimal('0.001'), rounding=ROUND_HALF_UP))\n",
    "    pval = tempDF2['AdjPval'].loc[xvar+'-vs-'+yvar]\n",
    "    below_limit = 0#Initialize\n",
    "    if pval==1.0:\n",
    "        pval_text = '1.0'\n",
    "    else:\n",
    "        if pval==0.0:#Due to smaller than the float minimum\n",
    "            pval = sys.float_info.min\n",
    "            print('P-value was smaller than the float minimum:', pval)\n",
    "            below_limit = 1\n",
    "        pval_text = f'{Decimal(str(pval)):.3E}'#Take more digits because rounding is bad here\n",
    "        significand, exponent = pval_text.split(sep='E-')\n",
    "        significand = str(Decimal(significand).quantize(Decimal('0.1'), rounding=ROUND_HALF_UP))\n",
    "        if significand=='10.0':\n",
    "            significand = '1.0'\n",
    "            exponent = str(int(exponent)-1)\n",
    "        if int(exponent)>2:\n",
    "            pval_text = significand+r'$\\times$'+'10'+r'$^{{-{0}}}$'.format(exponent)##Font is different in r'$ $'...\n",
    "        elif int(exponent)>0:\n",
    "            pval_text = '0.'+'0'*(int(exponent)-1)+significand.replace('.', '')\n",
    "        else:\n",
    "            pval_text = significand\n",
    "    if below_limit==1:\n",
    "        text = 'Pearson\\'s '+r'$r$'+' = '+r_text+'\\n'+r'$P$'+' < '+pval_text\n",
    "    else:\n",
    "        text = 'Pearson\\'s '+r'$r$'+' = '+r_text+'\\n'+r'$P$'+' = '+pval_text\n",
    "    p.axes[i, j].annotate(text, xy=(0.025, 1.0), xycoords='axes fraction',\n",
    "                          horizontalalignment='left', verticalalignment='top',\n",
    "                          multialignment='left', fontsize='x-small', color='k')\n",
    "pl = plt.legend(bbox_to_anchor=(0.8, 4.1), loc='upper right', title='BMI class')\n",
    "##Add sample size in lagend\n",
    "for row_i in range(len(pl.get_texts())):\n",
    "    bmi_class = pl.get_texts()[row_i].get_text()\n",
    "    count = len(tempDF.loc[tempDF['BaseBMI_class']==bmi_class])\n",
    "    pl.get_texts()[row_i].set_text(bmi_class+' ('+r'$n$'+' = '+f'{count:,}'+')')\n",
    "##Add xy label annotation\n",
    "label = r'$\\Delta$'+'BMI (predicted '+r'$-$'+' measured) [% BMI]'\n",
    "p.fig.text(x=0.545, y=0.0,#Manual adjustment\n",
    "           s=label, fontsize='large',\n",
    "           verticalalignment='top', horizontalalignment='center')\n",
    "p.fig.text(x=-0.01, y=0.545,#Manual adjustment\n",
    "           s=label, fontsize='large',\n",
    "           verticalalignment='center', horizontalalignment='right', rotation='vertical')\n",
    "##Save\n",
    "fileDir = './ExportFigures/'\n",
    "ipynbName = '220803_Multiomics-BMI-NatMed1stRevision_DeltaBMI-misclassification_'\n",
    "fileName = 'bBMI-BMI-difference.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()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Metabolomics space"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2-1. Standardization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Import the cleaned baseline omics dataframes\n",
    "fileDir = '../210104_Biological-BMI-paper/ExportData/'\n",
    "ipynbName = '210104_Biological-BMI-paper_RF-imputation_'\n",
    "fileName = 'baseline-metDF-with-RF-imputation.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",
    "##Drop BMI and covariates\n",
    "tempL = ['log_BaseBMI', 'Sex', 'BaseAge', 'PC1', 'PC2', 'PC3', 'PC4', 'PC5', 'Race']\n",
    "tempDF = tempDF.drop(columns=tempL)\n",
    "\n",
    "#Z-score transformation\n",
    "scaler = StandardScaler(copy=True, with_mean=True, with_std=True)\n",
    "tempA = scaler.fit_transform(tempDF)#Column direction\n",
    "tempDF = pd.DataFrame(data=tempA, index=tempDF.index, columns=tempDF.columns)\n",
    "\n",
    "display(tempDF.describe(include='all'))\n",
    "\n",
    "zscoreDF = tempDF"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2-2. Misclassification label"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Misclassification label dataframe\n",
    "labelDF = bmiDF.loc[bmiDF['BaseMetBMI_class']!='NotCalculated'][['BaseBMI_class', 'BaseMetBMI_class']]\n",
    "print('Sample size of Z-score DF:', len(zscoreDF))\n",
    "print('Sample size of label DF:', len(labelDF))\n",
    "\n",
    "#Misclassification\n",
    "tempL = []\n",
    "for row_n in labelDF.index.tolist():\n",
    "    if labelDF.loc[row_n, 'BaseBMI_class'] == labelDF.loc[row_n, 'BaseMetBMI_class']:\n",
    "        tempL.append('Matched')\n",
    "    else:\n",
    "        tempL.append('Mismatched')\n",
    "labelDF['Misclassification'] = tempL\n",
    "\n",
    "print('\\nSummary')\n",
    "tempL = ['Underweight', 'Normal', 'Overweight', 'Obese']\n",
    "for bmi_class in tempL:\n",
    "    print(' • BMI class: '+bmi_class)\n",
    "    tempDF = labelDF.loc[labelDF['BaseBMI_class']==bmi_class]\n",
    "    print('    - Sample size:', len(tempDF))\n",
    "    nMismatched = len(tempDF.loc[tempDF['Misclassification']=='Mismatched'])\n",
    "    print('    - Mismatched:', nMismatched, '(', nMismatched/len(tempDF)*100, '%)')\n",
    "    tempS = tempDF['BaseMetBMI_class'].value_counts()\n",
    "    tempDF = pd.DataFrame({'Count': tempS, 'Percentage':tempS/len(tempDF)*100})\n",
    "    print('    - Detail of MetBMI class:')\n",
    "    display(tempDF)\n",
    "    print('')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2-3. Color label for variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Prepare column color label based on the contribution in LASSO models\n",
    "fileDir = './ExportData/'\n",
    "ipynbName = '220802_Multiomics-BMI-NatMed1stRevision_BMI-LASSO-bcoef_'\n",
    "fileName = 'MetBMI-BothSex-OLS.tsv'\n",
    "tempDF = pd.read_csv(fileDir+ipynbName+fileName, sep='\\t').set_index('Variable')\n",
    "tempL = []\n",
    "for row_n in tempDF.index.tolist():\n",
    "    if row_n=='MetBMI model':\n",
    "        tempL.append('NaN')\n",
    "    else:\n",
    "        if tempDF.loc[row_n, 'LASSObcoef'] > 0:\n",
    "            tempL.append('Positive')\n",
    "        elif tempDF.loc[row_n, 'LASSObcoef'] < 0:\n",
    "            tempL.append('Negative')\n",
    "        else:\n",
    "            tempL.append('Error')\n",
    "tempDF['Contribution'] = tempL\n",
    "tempDF = tempDF.loc[tempDF.index!='MetBMI model']\n",
    "tempD = {'Positive':'tab:red', 'Negative':'tab:blue'}\n",
    "colorDF = pd.DataFrame({'Contribution\\nin LASSO':tempDF['Contribution'].map(tempD)})\n",
    "display(colorDF)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2-4. Hierarchical clustering: variables retained in ≥6 models and exhibiting ≥14% explained variance in BMI (top 15)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Import the cleaned OLS regression result for LASSO variable\n",
    "fileDir = './ExportData/'\n",
    "ipynbName = '220802_Multiomics-BMI-NatMed1stRevision_BMI-LASSO-bcoef_'\n",
    "fileName = 'MetBMI-BothSex-OLS.tsv'\n",
    "tempDF = pd.read_csv(fileDir+ipynbName+fileName, sep='\\t').set_index('Variable')\n",
    "tempDF = tempDF.loc[tempDF['LASSOnZeros']<=4]\n",
    "tempDF = tempDF.loc[tempDF['R2']>=14]\n",
    "#tempDF = tempDF.iloc[0:15]\n",
    "print('nrows:', len(tempDF))\n",
    "print('Adjusted p-value < 0.05:', len(tempDF.loc[tempDF['AdjPval']<0.05]))\n",
    "display(tempDF)\n",
    "\n",
    "#Extract the variables\n",
    "zscoreDF_select = zscoreDF.loc[:, zscoreDF.columns.isin(tempDF.index)]\n",
    "print('Confirm Z-score DF:', zscoreDF_select.shape)\n",
    "\n",
    "#Update color label\n",
    "var_colorDF = colorDF.loc[colorDF.index.isin(zscoreDF_select.columns)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('All:')\n",
    "#Visualization\n",
    "tempD1 = {'Underweight':'tab:blue', 'Normal':'tab:green', 'Overweight':'tab:orange', 'Obese':'tab:red'}\n",
    "tempD2 = {'Matched':'white', 'Mismatched':'black'}\n",
    "tempDF1 = pd.DataFrame({'BMI class':labelDF['BaseBMI_class'].map(tempD1),\n",
    "                        'MetBMI class':labelDF['BaseMetBMI_class'].map(tempD1),\n",
    "                        'Misclassification':labelDF['Misclassification'].map(tempD2)})\n",
    "t_start = time.time()\n",
    "sns.set(style='ticks', font='Arial', context='notebook')\n",
    "cm = sns.clustermap(zscoreDF_select, method='ward', metric='euclidean',\n",
    "                    row_cluster=True, col_cluster=True, row_linkage=None, col_linkage=None,\n",
    "                    row_colors=tempDF1, col_colors=var_colorDF,\n",
    "                    cbar_pos=(0.12, 0.01, 0.02, 0.1),\n",
    "                    figsize=(10, 30), **{'center':0, 'vmin':-2.58, 'vmax':2.58})\n",
    "cm.cax.set_title(r'$Z$'+'-score')\n",
    "hm = cm.ax_heatmap.get_position()\n",
    "rd = cm.ax_row_dendrogram.get_position()\n",
    "cd = cm.ax_col_dendrogram.get_position()\n",
    "cm.ax_heatmap.set_position([hm.x0, hm.y0, hm.width, hm.height])\n",
    "cm.ax_row_dendrogram.set_position([rd.x0+rd.width*0.5, rd.y0, rd.width*0.5, rd.height])\n",
    "cm.ax_col_dendrogram.set_position([cd.x0, cd.y0, cd.width, cd.height*0.25])\n",
    "plt.show()\n",
    "t_elapsed = time.time() - t_start\n",
    "print('• Elapsed time:', round(t_elapsed//60), 'min', round(t_elapsed%60, 1), 'sec')\n",
    "#Statistical test\n",
    "for nclusters in range(2, 6):\n",
    "    print('• nCluster: '+str(nclusters))\n",
    "    ##Perform hierarchical clustering for rows: public_client_id\n",
    "    tempDF1 = pd.DataFrame(hierarchy.linkage(zscoreDF_select, method='ward', metric='euclidean'))\n",
    "    ##Extract and merge cluster label\n",
    "    tempS = pd.Series(hierarchy.fcluster(tempDF1, nclusters, criterion='maxclust'),\n",
    "                      index=zscoreDF_select.index, name='Cluster', dtype=str)\n",
    "    tempS = 'Cluster_'+tempS\n",
    "    tempDF = pd.merge(labelDF, tempS, left_index=True, right_index=True)\n",
    "    ##To check which cluster corresponds to sns.clustermap\n",
    "    tempL1 = []\n",
    "    for row_i in range(len(cm.ax_heatmap.get_yticklabels())):\n",
    "        tempL1.append(cm.ax_heatmap.get_yticklabels()[row_i].get_text())\n",
    "    tempDF1 = tempDF.loc[tempL1]\n",
    "    print('  - For checking cluster label')\n",
    "    for cluster in tempDF['Cluster'].unique().tolist():\n",
    "        print('     - '+cluster+' includes', tempDF1.loc[tempDF1['Cluster']==cluster].index.tolist()[0:5])\n",
    "    ##Summary\n",
    "    tempDF = pd.crosstab(tempDF['Misclassification'], tempDF['Cluster'], margins=True)\n",
    "    print('  - Misclassification vs. Cluster')\n",
    "    display(tempDF)\n",
    "    ##Fisher's exact test for each cluster\n",
    "    tempDF = tempDF.drop(index=['All'])\n",
    "    tempS = pd.Series(index=tempDF.drop(columns=['All']).columns, name='Pval')#For p-values\n",
    "    for cluster in tempS.index.tolist():\n",
    "        #Prepare 2x2 matrix\n",
    "        tempDF1 = pd.DataFrame({cluster:tempDF[cluster],\n",
    "                                'Others':tempDF['All']-tempDF[cluster]})\n",
    "        if tempDF1.shape == (2, 2):\n",
    "            #Fisher's exact test\n",
    "            tempS[cluster] = stats.fisher_exact(tempDF1, alternative='two-sided')[1]\n",
    "        else:\n",
    "            tempS[cluster] = np.nan\n",
    "    ##Multiple tests correction\n",
    "    tempDF = pd.DataFrame({'Pval':tempS,\n",
    "                           'AdjPval':multi.multipletests(tempS, alpha=0.05, method='bonferroni',\n",
    "                                                         is_sorted=False, returnsorted=False)[1]})\n",
    "    print('  - Fisher\\'s exact test (two-sided, Bonferroni correction): cluster vs. the other clusters')\n",
    "    display(tempDF)\n",
    "    print('')\n",
    "print('')\n",
    "\n",
    "#BMI class\n",
    "tempL = ['Normal', 'Overweight', 'Obese']\n",
    "for bmi_class in tempL:\n",
    "    print('BMI class: '+bmi_class)\n",
    "    labelDF_temp = labelDF.loc[labelDF['BaseBMI_class']==bmi_class]\n",
    "    tempDF = zscoreDF_select.loc[zscoreDF_select.index.isin(labelDF_temp.index)]\n",
    "    #Visualization\n",
    "    tempD1 = {'Underweight':'tab:blue', 'Normal':'tab:green', 'Overweight':'tab:orange', 'Obese':'tab:red'}\n",
    "    tempD2 = {'Matched':'white', 'Mismatched':'black'}\n",
    "    tempDF1 = pd.DataFrame({'BMI class':labelDF_temp['BaseBMI_class'].map(tempD1),\n",
    "                            'MetBMI class':labelDF_temp['BaseMetBMI_class'].map(tempD1),\n",
    "                            'Misclassification':labelDF_temp['Misclassification'].map(tempD2)})\n",
    "    t_start = time.time()\n",
    "    sns.set(style='ticks', font='Arial', context='notebook')\n",
    "    cm = sns.clustermap(tempDF, method='ward', metric='euclidean',\n",
    "                        row_cluster=True, col_cluster=True, row_linkage=None, col_linkage=None,\n",
    "                        row_colors=tempDF1, col_colors=var_colorDF,\n",
    "                        cbar_pos=(0.12, 0.02, 0.02, 0.2),\n",
    "                        figsize=(10, 15), **{'center':0, 'vmin':-2.58, 'vmax':2.58})\n",
    "    cm.cax.set_title(r'$Z$'+'-score')\n",
    "    hm = cm.ax_heatmap.get_position()\n",
    "    rd = cm.ax_row_dendrogram.get_position()\n",
    "    cd = cm.ax_col_dendrogram.get_position()\n",
    "    cm.ax_heatmap.set_position([hm.x0, hm.y0, hm.width, hm.height])\n",
    "    cm.ax_row_dendrogram.set_position([rd.x0+rd.width*0.5, rd.y0, rd.width*0.5, rd.height])\n",
    "    cm.ax_col_dendrogram.set_position([cd.x0, cd.y0, cd.width, cd.height*0.25])\n",
    "    plt.show()\n",
    "    t_elapsed = time.time() - t_start\n",
    "    print('• Elapsed time:', round(t_elapsed//60), 'min', round(t_elapsed%60, 1), 'sec')\n",
    "    #Statistical test\n",
    "    for nclusters in range(2, 6):\n",
    "        print('• nCluster: '+str(nclusters))\n",
    "        ##Perform hierarchical clustering for rows: public_client_id\n",
    "        tempDF1 = pd.DataFrame(hierarchy.linkage(tempDF, method='ward', metric='euclidean'))\n",
    "        ##Extract and merge cluster label\n",
    "        tempS = pd.Series(hierarchy.fcluster(tempDF1, nclusters, criterion='maxclust'),\n",
    "                          index=tempDF.index, name='Cluster', dtype=str)\n",
    "        tempS = 'Cluster_'+tempS\n",
    "        tempDF2 = pd.merge(labelDF_temp, tempS, left_index=True, right_index=True)\n",
    "        ##To check which cluster corresponds to sns.clustermap\n",
    "        tempL1 = []\n",
    "        for row_i in range(len(cm.ax_heatmap.get_yticklabels())):\n",
    "            tempL1.append(cm.ax_heatmap.get_yticklabels()[row_i].get_text())\n",
    "        tempDF1 = tempDF2.loc[tempL1]\n",
    "        print('  - For checking cluster label')\n",
    "        for cluster in tempDF2['Cluster'].unique().tolist():\n",
    "            print('     - '+cluster+' includes', tempDF1.loc[tempDF1['Cluster']==cluster].index.tolist()[0:5])\n",
    "        ##Summary\n",
    "        tempDF1 = pd.crosstab(tempDF2['Misclassification'], tempDF2['Cluster'], margins=True)\n",
    "        print('  - Misclassification vs. Cluster')\n",
    "        display(tempDF1)\n",
    "        ##Fisher's exact test for each cluster\n",
    "        tempDF1 = tempDF1.drop(index=['All'])\n",
    "        tempS = pd.Series(index=tempDF1.drop(columns=['All']).columns, name='Pval')#For p-values\n",
    "        for cluster in tempS.index.tolist():\n",
    "            #Prepare 2x2 matrix\n",
    "            tempDF2 = pd.DataFrame({cluster:tempDF1[cluster],\n",
    "                                    'Others':tempDF1['All']-tempDF1[cluster]})\n",
    "            if tempDF2.shape == (2, 2):\n",
    "                #Fisher's exact test\n",
    "                tempS[cluster] = stats.fisher_exact(tempDF2, alternative='two-sided')[1]\n",
    "            else:\n",
    "                tempS[cluster] = np.nan\n",
    "        ##Multiple tests correction\n",
    "        tempDF1 = pd.DataFrame({'Pval':tempS,\n",
    "                               'AdjPval':multi.multipletests(tempS, alpha=0.05, method='bonferroni',\n",
    "                                                             is_sorted=False, returnsorted=False)[1]})\n",
    "        print('  - Fisher\\'s exact test (two-sided, Bonferroni correction): cluster vs. the other clusters')\n",
    "        display(tempDF1)\n",
    "        print('')\n",
    "    print('')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Plot for paper\n",
    "tempL = ['Normal', 'Obese']\n",
    "for bmi_class in tempL:\n",
    "    print('BMI class: '+bmi_class)\n",
    "    labelDF_temp = labelDF.loc[labelDF['BaseBMI_class']==bmi_class]\n",
    "    tempDF = zscoreDF_select.loc[zscoreDF_select.index.isin(labelDF_temp.index)]\n",
    "    #Visualization\n",
    "    tempD1 = {'Underweight':'tab:blue', 'Normal':'tab:green', 'Overweight':'tab:orange', 'Obese':'tab:red'}\n",
    "    tempD2 = {'Matched':'white', 'Mismatched':'black'}\n",
    "    tempDF1 = pd.DataFrame({'BMI class':labelDF_temp['BaseBMI_class'].map(tempD1),\n",
    "                            'MetBMI class':labelDF_temp['BaseMetBMI_class'].map(tempD1),\n",
    "                            'Misclassification':labelDF_temp['Misclassification'].map(tempD2)})\n",
    "    t_start = time.time()\n",
    "    sns.set(style='ticks', font='Arial', context='talk')\n",
    "    cm = sns.clustermap(tempDF, method='ward', metric='euclidean',\n",
    "                        row_cluster=True, col_cluster=True, row_linkage=None, col_linkage=None,\n",
    "                        row_colors=tempDF1, col_colors=var_colorDF, xticklabels=True, yticklabels=False,\n",
    "                        dendrogram_ratio=(0.2, 0.1), colors_ratio=(0.045, 0.03),\n",
    "                        cbar_pos=(0.85, 0.43, 0.3, 0.02), cbar_kws={\"orientation\": \"horizontal\"},\n",
    "                        figsize=(8, 13.45), **{'center':0, 'vmin':-2.58, 'vmax':2.58})\n",
    "    cm.cax.set_title(r'$Z$'+'-score', size='medium', verticalalignment='bottom')\n",
    "    cm.cax.tick_params(labelsize='small')\n",
    "    bottom, top = cm.ax_heatmap.get_ylim()\n",
    "    cm.ax_heatmap.set_ylim(bottom + 0.5, top - 0.5)##To avoid half cut of first and last rows\n",
    "    hm = cm.ax_heatmap.get_position()\n",
    "    rd = cm.ax_row_dendrogram.get_position()\n",
    "    cd = cm.ax_col_dendrogram.get_position()\n",
    "    cm.ax_heatmap.set_position([hm.x0, hm.y0, hm.width, hm.height])\n",
    "    cm.ax_row_dendrogram.set_position([rd.x0+rd.width*0, rd.y0, rd.width, rd.height])\n",
    "    cm.ax_col_dendrogram.set_position([cd.x0, cd.y0, cd.width, cd.height*0.5])\n",
    "    cm.ax_heatmap.set_xlabel('')\n",
    "    #cm.ax_heatmap.set_ylabel('Participant ('+r'$n$ = '+f'{len(tempDF):,}'+')')\n",
    "    cm.ax_heatmap.set_ylabel('')\n",
    "    #cm.ax_col_colors.yaxis.set_ticks_position('left')\n",
    "    ##row/column color bar legend (axis is same with cm.cax!)\n",
    "    row_legend1 = mpatches.Patch(color='tab:green', label='Normal')\n",
    "    row_legend2 = mpatches.Patch(color='tab:orange', label='Overweight')\n",
    "    row_legend3 = mpatches.Patch(color='tab:red', label='Obese')\n",
    "    legend1 = plt.legend(handles=[row_legend1, row_legend2, row_legend3], fontsize='medium',\n",
    "                         title='BMI/MetBMI class', title_fontsize='medium',\n",
    "                         bbox_to_anchor=(0.5, 0), loc='upper center', borderaxespad=3.5, frameon=False)\n",
    "    plt.gca().add_artist(legend1)\n",
    "    row_legend = mpatches.Patch(color='black', label='Mismatched')\n",
    "    legend2 = plt.legend(handles=[row_legend], fontsize='medium',\n",
    "                         title='Misclassification', title_fontsize='medium',\n",
    "                         bbox_to_anchor=(0.5, 0), loc='upper center', borderaxespad=9.5, frameon=False)\n",
    "    plt.gca().add_artist(legend2)\n",
    "    col_legend1 = mpatches.Patch(color='tab:red', label='Positive')\n",
    "    col_legend2 = mpatches.Patch(color='tab:blue', label='Negative')\n",
    "    plt.legend(handles=[col_legend1, col_legend2], fontsize='medium',\n",
    "               title='Contribution in LASSO', title_fontsize='medium',\n",
    "               bbox_to_anchor=(0.5, 0), loc='upper center', borderaxespad=15, frameon=False)\n",
    "    ##Save\n",
    "    fileDir = './ExportFigures/'\n",
    "    ipynbName = '220803_Multiomics-BMI-NatMed1stRevision_DeltaBMI-misclassification_'\n",
    "    fileName = 'clusters-15metabolites-'+bmi_class+'.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",
    "    t_elapsed = time.time() - t_start\n",
    "    print('• Elapsed time:', round(t_elapsed//60), 'min', round(t_elapsed%60, 1), 'sec')\n",
    "    print('')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Proteomics space"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3-1. Standardization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Import the cleaned baseline omics dataframes\n",
    "fileDir = '../210104_Biological-BMI-paper/ExportData/'\n",
    "ipynbName = '210104_Biological-BMI-paper_RF-imputation_'\n",
    "fileName = 'baseline-protDF-with-RF-imputation.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",
    "##Drop BMI and covariates\n",
    "tempL = ['log_BaseBMI', 'Sex', 'BaseAge', 'PC1', 'PC2', 'PC3', 'PC4', 'PC5', 'Race']\n",
    "tempDF = tempDF.drop(columns=tempL)\n",
    "\n",
    "#Z-score transformation\n",
    "scaler = StandardScaler(copy=True, with_mean=True, with_std=True)\n",
    "tempA = scaler.fit_transform(tempDF)#Column direction\n",
    "tempDF = pd.DataFrame(data=tempA, index=tempDF.index, columns=tempDF.columns)\n",
    "\n",
    "display(tempDF.describe(include='all'))\n",
    "\n",
    "zscoreDF = tempDF"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3-2. Misclassification label"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Misclassification label dataframe\n",
    "labelDF = bmiDF.loc[bmiDF['BaseProtBMI_class']!='NotCalculated'][['BaseBMI_class', 'BaseProtBMI_class']]\n",
    "print('Sample size of Z-score DF:', len(zscoreDF))\n",
    "print('Sample size of label DF:', len(labelDF))\n",
    "\n",
    "#Misclassification\n",
    "tempL = []\n",
    "for row_n in labelDF.index.tolist():\n",
    "    if labelDF.loc[row_n, 'BaseBMI_class'] == labelDF.loc[row_n, 'BaseProtBMI_class']:\n",
    "        tempL.append('Matched')\n",
    "    else:\n",
    "        tempL.append('Mismatched')\n",
    "labelDF['Misclassification'] = tempL\n",
    "\n",
    "print('\\nSummary')\n",
    "tempL = ['Underweight', 'Normal', 'Overweight', 'Obese']\n",
    "for bmi_class in tempL:\n",
    "    print(' • BMI class: '+bmi_class)\n",
    "    tempDF = labelDF.loc[labelDF['BaseBMI_class']==bmi_class]\n",
    "    print('    - Sample size:', len(tempDF))\n",
    "    nMismatched = len(tempDF.loc[tempDF['Misclassification']=='Mismatched'])\n",
    "    print('    - Mismatched:', nMismatched, '(', nMismatched/len(tempDF)*100, '%)')\n",
    "    tempS = tempDF['BaseProtBMI_class'].value_counts()\n",
    "    tempDF = pd.DataFrame({'Count': tempS, 'Percentage':tempS/len(tempDF)*100})\n",
    "    print('    - Detail of ProtBMI class:')\n",
    "    display(tempDF)\n",
    "    print('')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3-3. Color label for variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Prepare column color label based on the contribution in LASSO models\n",
    "fileDir = './ExportData/'\n",
    "ipynbName = '220802_Multiomics-BMI-NatMed1stRevision_BMI-LASSO-bcoef_'\n",
    "fileName = 'ProtBMI-BothSex-OLS.tsv'\n",
    "tempDF = pd.read_csv(fileDir+ipynbName+fileName, sep='\\t').set_index('Variable')\n",
    "tempL = []\n",
    "for row_n in tempDF.index.tolist():\n",
    "    if row_n=='ProtBMI model':\n",
    "        tempL.append('NaN')\n",
    "    else:\n",
    "        if tempDF.loc[row_n, 'LASSObcoef'] > 0:\n",
    "            tempL.append('Positive')\n",
    "        elif tempDF.loc[row_n, 'LASSObcoef'] < 0:\n",
    "            tempL.append('Negative')\n",
    "        else:\n",
    "            tempL.append('Error')\n",
    "tempDF['Contribution'] = tempL\n",
    "tempDF = tempDF.loc[tempDF.index!='ProtBMI model']\n",
    "tempD = {'Positive':'tab:red', 'Negative':'tab:blue'}\n",
    "colorDF = pd.DataFrame({'Contribution\\nin LASSO':tempDF['Contribution'].map(tempD)})\n",
    "display(colorDF)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3-4. Hierarchical clustering: variables retained in ≥8 models and exhibiting ≥10% explained variance in BMI (top 15)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Import the cleaned OLS regression result for LASSO variable\n",
    "fileDir = './ExportData/'\n",
    "ipynbName = '220802_Multiomics-BMI-NatMed1stRevision_BMI-LASSO-bcoef_'\n",
    "fileName = 'ProtBMI-BothSex-OLS.tsv'\n",
    "tempDF = pd.read_csv(fileDir+ipynbName+fileName, sep='\\t').set_index('Variable')\n",
    "tempDF = tempDF.loc[tempDF['LASSOnZeros']<=2]\n",
    "tempDF = tempDF.loc[tempDF['R2']>=10]\n",
    "#tempDF = tempDF.iloc[0:15]\n",
    "print('nrows:', len(tempDF))\n",
    "print('Adjusted p-value < 0.05:', len(tempDF.loc[tempDF['AdjPval']<0.05]))\n",
    "display(tempDF)\n",
    "\n",
    "#Extract the variables\n",
    "zscoreDF_select = zscoreDF.loc[:, zscoreDF.columns.isin(tempDF.index)]\n",
    "print('Confirm Z-score DF:', zscoreDF_select.shape)\n",
    "\n",
    "#Update color label\n",
    "var_colorDF = colorDF.loc[colorDF.index.isin(zscoreDF_select.columns)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('All:')\n",
    "#Visualization\n",
    "tempD1 = {'Underweight':'tab:blue', 'Normal':'tab:green', 'Overweight':'tab:orange', 'Obese':'tab:red'}\n",
    "tempD2 = {'Matched':'white', 'Mismatched':'black'}\n",
    "tempDF1 = pd.DataFrame({'BMI class':labelDF['BaseBMI_class'].map(tempD1),\n",
    "                        'ProtBMI class':labelDF['BaseProtBMI_class'].map(tempD1),\n",
    "                        'Misclassification':labelDF['Misclassification'].map(tempD2)})\n",
    "t_start = time.time()\n",
    "sns.set(style='ticks', font='Arial', context='notebook')\n",
    "cm = sns.clustermap(zscoreDF_select, method='ward', metric='euclidean',\n",
    "                    row_cluster=True, col_cluster=True, row_linkage=None, col_linkage=None,\n",
    "                    row_colors=tempDF1, col_colors=var_colorDF,\n",
    "                    cbar_pos=(0.12, 0.825, 0.15, 0.01), cbar_kws={\"orientation\": \"horizontal\"},\n",
    "                    figsize=(10, 24), **{'center':0, 'vmin':-2.58, 'vmax':2.58})\n",
    "cm.cax.set_title(r'$Z$'+'-score')\n",
    "hm = cm.ax_heatmap.get_position()\n",
    "rd = cm.ax_row_dendrogram.get_position()\n",
    "cd = cm.ax_col_dendrogram.get_position()\n",
    "cm.ax_heatmap.set_position([hm.x0, hm.y0, hm.width, hm.height])\n",
    "cm.ax_row_dendrogram.set_position([rd.x0+rd.width*0.5, rd.y0, rd.width*0.5, rd.height])\n",
    "cm.ax_col_dendrogram.set_position([cd.x0, cd.y0, cd.width, cd.height*0.25])\n",
    "plt.show()\n",
    "t_elapsed = time.time() - t_start\n",
    "print('• Elapsed time:', round(t_elapsed//60), 'min', round(t_elapsed%60, 1), 'sec')\n",
    "#Statistical test\n",
    "for nclusters in range(2, 6):\n",
    "    print('• nCluster: '+str(nclusters))\n",
    "    ##Perform hierarchical clustering for rows: public_client_id\n",
    "    tempDF1 = pd.DataFrame(hierarchy.linkage(zscoreDF_select, method='ward', metric='euclidean'))\n",
    "    ##Extract and merge cluster label\n",
    "    tempS = pd.Series(hierarchy.fcluster(tempDF1, nclusters, criterion='maxclust'),\n",
    "                      index=zscoreDF_select.index, name='Cluster', dtype=str)\n",
    "    tempS = 'Cluster_'+tempS\n",
    "    tempDF = pd.merge(labelDF, tempS, left_index=True, right_index=True)\n",
    "    ##To check which cluster corresponds to sns.clustermap\n",
    "    tempL1 = []\n",
    "    for row_i in range(len(cm.ax_heatmap.get_yticklabels())):\n",
    "        tempL1.append(cm.ax_heatmap.get_yticklabels()[row_i].get_text())\n",
    "    tempDF1 = tempDF.loc[tempL1]\n",
    "    print('  - For checking cluster label')\n",
    "    for cluster in tempDF['Cluster'].unique().tolist():\n",
    "        print('     - '+cluster+' includes', tempDF1.loc[tempDF1['Cluster']==cluster].index.tolist()[0:5])\n",
    "    ##Summary\n",
    "    tempDF = pd.crosstab(tempDF['Misclassification'], tempDF['Cluster'], margins=True)\n",
    "    print('  - Misclassification vs. Cluster')\n",
    "    display(tempDF)\n",
    "    ##Fisher's exact test for each cluster\n",
    "    tempDF = tempDF.drop(index=['All'])\n",
    "    tempS = pd.Series(index=tempDF.drop(columns=['All']).columns, name='Pval')#For p-values\n",
    "    for cluster in tempS.index.tolist():\n",
    "        #Prepare 2x2 matrix\n",
    "        tempDF1 = pd.DataFrame({cluster:tempDF[cluster],\n",
    "                                'Others':tempDF['All']-tempDF[cluster]})\n",
    "        if tempDF1.shape == (2, 2):\n",
    "            #Fisher's exact test\n",
    "            tempS[cluster] = stats.fisher_exact(tempDF1, alternative='two-sided')[1]\n",
    "        else:\n",
    "            tempS[cluster] = np.nan\n",
    "    ##Multiple tests correction\n",
    "    tempDF = pd.DataFrame({'Pval':tempS,\n",
    "                           'AdjPval':multi.multipletests(tempS, alpha=0.05, method='bonferroni',\n",
    "                                                         is_sorted=False, returnsorted=False)[1]})\n",
    "    print('  - Fisher\\'s exact test (two-sided, Bonferroni correction): cluster vs. the other clusters')\n",
    "    display(tempDF)\n",
    "    print('')\n",
    "print('')\n",
    "\n",
    "#BMI class\n",
    "tempL = ['Normal', 'Overweight', 'Obese']\n",
    "for bmi_class in tempL:\n",
    "    print('BMI class: '+bmi_class)\n",
    "    labelDF_temp = labelDF.loc[labelDF['BaseBMI_class']==bmi_class]\n",
    "    tempDF = zscoreDF_select.loc[zscoreDF_select.index.isin(labelDF_temp.index)]\n",
    "    #Visualization\n",
    "    tempD1 = {'Underweight':'tab:blue', 'Normal':'tab:green', 'Overweight':'tab:orange', 'Obese':'tab:red'}\n",
    "    tempD2 = {'Matched':'white', 'Mismatched':'black'}\n",
    "    tempDF1 = pd.DataFrame({'BMI class':labelDF_temp['BaseBMI_class'].map(tempD1),\n",
    "                            'ProtBMI class':labelDF_temp['BaseProtBMI_class'].map(tempD1),\n",
    "                            'Misclassification':labelDF_temp['Misclassification'].map(tempD2)})\n",
    "    t_start = time.time()\n",
    "    sns.set(style='ticks', font='Arial', context='notebook')\n",
    "    cm = sns.clustermap(tempDF, method='ward', metric='euclidean',\n",
    "                        row_cluster=True, col_cluster=True, row_linkage=None, col_linkage=None,\n",
    "                        row_colors=tempDF1, col_colors=var_colorDF,\n",
    "                        cbar_pos=(0.12, 0.825, 0.15, 0.02), cbar_kws={\"orientation\": \"horizontal\"},\n",
    "                        figsize=(10, 12), **{'center':0, 'vmin':-2.58, 'vmax':2.58})\n",
    "    cm.cax.set_title(r'$Z$'+'-score')\n",
    "    hm = cm.ax_heatmap.get_position()\n",
    "    rd = cm.ax_row_dendrogram.get_position()\n",
    "    cd = cm.ax_col_dendrogram.get_position()\n",
    "    cm.ax_heatmap.set_position([hm.x0, hm.y0, hm.width, hm.height])\n",
    "    cm.ax_row_dendrogram.set_position([rd.x0+rd.width*0.5, rd.y0, rd.width*0.5, rd.height])\n",
    "    cm.ax_col_dendrogram.set_position([cd.x0, cd.y0, cd.width, cd.height*0.25])\n",
    "    plt.show()\n",
    "    t_elapsed = time.time() - t_start\n",
    "    print('• Elapsed time:', round(t_elapsed//60), 'min', round(t_elapsed%60, 1), 'sec')\n",
    "    #Statistical test\n",
    "    for nclusters in range(2, 6):\n",
    "        print('• nCluster: '+str(nclusters))\n",
    "        ##Perform hierarchical clustering for rows: public_client_id\n",
    "        tempDF1 = pd.DataFrame(hierarchy.linkage(tempDF, method='ward', metric='euclidean'))\n",
    "        ##Extract and merge cluster label\n",
    "        tempS = pd.Series(hierarchy.fcluster(tempDF1, nclusters, criterion='maxclust'),\n",
    "                          index=tempDF.index, name='Cluster', dtype=str)\n",
    "        tempS = 'Cluster_'+tempS\n",
    "        tempDF2 = pd.merge(labelDF_temp, tempS, left_index=True, right_index=True)\n",
    "        ##To check which cluster corresponds to sns.clustermap\n",
    "        tempL1 = []\n",
    "        for row_i in range(len(cm.ax_heatmap.get_yticklabels())):\n",
    "            tempL1.append(cm.ax_heatmap.get_yticklabels()[row_i].get_text())\n",
    "        tempDF1 = tempDF2.loc[tempL1]\n",
    "        print('  - For checking cluster label')\n",
    "        for cluster in tempDF2['Cluster'].unique().tolist():\n",
    "            print('     - '+cluster+' includes', tempDF1.loc[tempDF1['Cluster']==cluster].index.tolist()[0:5])\n",
    "        ##Summary\n",
    "        tempDF1 = pd.crosstab(tempDF2['Misclassification'], tempDF2['Cluster'], margins=True)\n",
    "        print('  - Misclassification vs. Cluster')\n",
    "        display(tempDF1)\n",
    "        ##Fisher's exact test for each cluster\n",
    "        tempDF1 = tempDF1.drop(index=['All'])\n",
    "        tempS = pd.Series(index=tempDF1.drop(columns=['All']).columns, name='Pval')#For p-values\n",
    "        for cluster in tempS.index.tolist():\n",
    "            #Prepare 2x2 matrix\n",
    "            tempDF2 = pd.DataFrame({cluster:tempDF1[cluster],\n",
    "                                    'Others':tempDF1['All']-tempDF1[cluster]})\n",
    "            if tempDF2.shape == (2, 2):\n",
    "                #Fisher's exact test\n",
    "                tempS[cluster] = stats.fisher_exact(tempDF2, alternative='two-sided')[1]\n",
    "            else:\n",
    "                tempS[cluster] = np.nan\n",
    "        ##Multiple tests correction\n",
    "        tempDF1 = pd.DataFrame({'Pval':tempS,\n",
    "                               'AdjPval':multi.multipletests(tempS, alpha=0.05, method='bonferroni',\n",
    "                                                             is_sorted=False, returnsorted=False)[1]})\n",
    "        print('  - Fisher\\'s exact test (two-sided, Bonferroni correction): cluster vs. the other clusters')\n",
    "        display(tempDF1)\n",
    "        print('')\n",
    "    print('')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Plot for paper\n",
    "tempL = ['Normal', 'Obese']\n",
    "for bmi_class in tempL:\n",
    "    print('BMI class: '+bmi_class)\n",
    "    labelDF_temp = labelDF.loc[labelDF['BaseBMI_class']==bmi_class]\n",
    "    tempDF = zscoreDF_select.loc[zscoreDF_select.index.isin(labelDF_temp.index)]\n",
    "    #Visualization\n",
    "    tempD1 = {'Underweight':'tab:blue', 'Normal':'tab:green', 'Overweight':'tab:orange', 'Obese':'tab:red'}\n",
    "    tempD2 = {'Matched':'white', 'Mismatched':'black'}\n",
    "    tempDF1 = pd.DataFrame({'BMI class':labelDF_temp['BaseBMI_class'].map(tempD1),\n",
    "                            'ProtBMI class':labelDF_temp['BaseProtBMI_class'].map(tempD1),\n",
    "                            'Misclassification':labelDF_temp['Misclassification'].map(tempD2)})\n",
    "    t_start = time.time()\n",
    "    sns.set(style='ticks', font='Arial', context='talk')\n",
    "    cm = sns.clustermap(tempDF, method='ward', metric='euclidean',\n",
    "                        row_cluster=True, col_cluster=True, row_linkage=None, col_linkage=None,\n",
    "                        row_colors=tempDF1, col_colors=var_colorDF, xticklabels=True, yticklabels=False,\n",
    "                        dendrogram_ratio=(0.2, 0.1), colors_ratio=(0.045, 0.03),\n",
    "                        cbar_pos=(0.85, 0.575, 0.3, 0.02), cbar_kws={\"orientation\": \"horizontal\"},\n",
    "                        figsize=(8, 10), **{'center':0, 'vmin':-2.58, 'vmax':2.58})\n",
    "    cm.cax.set_title(r'$Z$'+'-score', size='medium', verticalalignment='bottom')\n",
    "    cm.cax.tick_params(labelsize='small')\n",
    "    bottom, top = cm.ax_heatmap.get_ylim()\n",
    "    cm.ax_heatmap.set_ylim(bottom + 0.5, top - 0.5)##To avoid half cut of first and last rows\n",
    "    hm = cm.ax_heatmap.get_position()\n",
    "    rd = cm.ax_row_dendrogram.get_position()\n",
    "    cd = cm.ax_col_dendrogram.get_position()\n",
    "    cm.ax_heatmap.set_position([hm.x0, hm.y0, hm.width, hm.height])\n",
    "    cm.ax_row_dendrogram.set_position([rd.x0+rd.width*0, rd.y0, rd.width, rd.height])\n",
    "    cm.ax_col_dendrogram.set_position([cd.x0, cd.y0, cd.width, cd.height*0.5])\n",
    "    cm.ax_heatmap.set_xlabel('')\n",
    "    #cm.ax_heatmap.set_ylabel('Participant ('+r'$n$ = '+f'{len(tempDF):,}'+')')\n",
    "    cm.ax_heatmap.set_ylabel('')\n",
    "    #cm.ax_col_colors.yaxis.set_ticks_position('left')\n",
    "    ##row/column color bar legend (axis is same with cm.cax!)\n",
    "    row_legend1 = mpatches.Patch(color='tab:green', label='Normal')\n",
    "    row_legend2 = mpatches.Patch(color='tab:orange', label='Overweight')\n",
    "    row_legend3 = mpatches.Patch(color='tab:red', label='Obese')\n",
    "    legend1 = plt.legend(handles=[row_legend1, row_legend2, row_legend3], fontsize='medium',\n",
    "                         title='BMI/ProtBMI class', title_fontsize='medium',\n",
    "                         bbox_to_anchor=(0.5, 0), loc='upper center', borderaxespad=3.5, frameon=False)\n",
    "    plt.gca().add_artist(legend1)\n",
    "    row_legend = mpatches.Patch(color='black', label='Mismatched')\n",
    "    legend2 = plt.legend(handles=[row_legend], fontsize='medium',\n",
    "                         title='Misclassification', title_fontsize='medium',\n",
    "                         bbox_to_anchor=(0.5, 0), loc='upper center', borderaxespad=9.5, frameon=False)\n",
    "    plt.gca().add_artist(legend2)\n",
    "    col_legend1 = mpatches.Patch(color='tab:red', label='Positive')\n",
    "    col_legend2 = mpatches.Patch(color='tab:blue', label='Negative')\n",
    "    plt.legend(handles=[col_legend1, col_legend2], fontsize='medium',\n",
    "               title='Contribution in LASSO', title_fontsize='medium',\n",
    "               bbox_to_anchor=(0.5, 0), loc='upper center', borderaxespad=15, frameon=False)\n",
    "    ##Save\n",
    "    fileDir = './ExportFigures/'\n",
    "    ipynbName = '220803_Multiomics-BMI-NatMed1stRevision_DeltaBMI-misclassification_'\n",
    "    fileName = 'clusters-15proteins-'+bmi_class+'.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",
    "    t_elapsed = time.time() - t_start\n",
    "    print('• Elapsed time:', round(t_elapsed//60), 'min', round(t_elapsed%60, 1), 'sec')\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"
  },
  "widgets": {
   "application/vnd.jupyter.widget-state+json": {
    "state": {},
    "version_major": 2,
    "version_minor": 0
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}