Code09_BMI-regression_numeric-features.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Multiomics BMI Paper — Regression Analysis of Measured and Biological BMIs on Numeric Features"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***by Kengo Watanabe***  \n",
    "\n",
    "This Jupyter Notebook (with Python 3 kernel) regressed the measured and omics-inferred BMIs independently on each of the available numeric physiological features (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 numeric physiological features (preprocessed): 210106_Biological-BMI-paper_data-cleaning_OLS-regression_numDF.tsv  \n",
    "* Numeric feature metadata: 210727_feature-label-correspondence.csv  \n",
    "\n",
    "Output figures and tables:  \n",
    "* Figure 1e  \n",
    "* Table for Supplementary Data 4  \n",
    "\n",
    "Original notebook (memo for my future tracing):  \n",
    "* dalek:\\[JupyterLab HOME\\]/220621_Multiomics-BMI-NatMedRevision/220803_Multiomics-BMI-NatMed1stRevision_BMI-regression_NumericFeatures.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 sklearn.preprocessing import StandardScaler\n",
    "import statsmodels.formula.api as smf\n",
    "from statsmodels.stats import multitest as multi\n",
    "\n",
    "!conda list"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Processing for OLS linear regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1-1. Prepare datasets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Import cleaned table for baseline measured and biological BMIs\n",
    "fileDir = './ExportData/'\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",
    "\n",
    "#Clean to handle easier in this notebook\n",
    "tempDF.columns = tempDF.columns.str.replace('Base', '')\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": [
    "#Import the cleaned DF of baseline numeric features\n",
    "fileDir = '../210104_Biological-BMI-paper/ExportData/'\n",
    "ipynbName = '210106_Biological-BMI-paper_data-cleaning_OLS-regression_'\n",
    "fileName = 'numDF.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",
    "#Clean to handle easier in this notebook\n",
    "tempDF = tempDF.rename(columns={'BaseAge':'Age'})\n",
    "tempDF = tempDF.drop(columns=['Race'])#Race has NaN in this cohort\n",
    "\n",
    "#Eliminate the exceptional feature that should have been eliminated during data cleaning\n",
    "feature = 'activities.calories'\n",
    "tempDF = tempDF.drop(columns=[feature])\n",
    "\n",
    "display(tempDF)\n",
    "print('NaN in DF:', tempDF.isnull().to_numpy().sum(axis=None))\n",
    "\n",
    "numDF = tempDF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Prepare the shortened feature label for visualization\n",
    "fileDir = '../210104_Biological-BMI-paper/ImportData/'\n",
    "fileName = '210727_feature-label-correspondence.csv'\n",
    "tempDF = pd.read_csv(fileDir+fileName)\n",
    "tempDF = tempDF.rename(columns={'OriginalName':'FeatureID', 'LabelForFigure':'FeatureLabel'})\n",
    "tempDF = tempDF.set_index('FeatureID')\n",
    "\n",
    "display(tempDF)\n",
    "print(' -> Unique label:', len(tempDF['FeatureLabel'].unique()))\n",
    "\n",
    "numDF_meta = tempDF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Define the list of covariates\n",
    "covarL = ['Sex', 'Age', 'PC1', 'PC2', 'PC3', 'PC4', 'PC5']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1-2. Split DF for each regression, and then drop NaN"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Sample size is different b/w numeric fetuares.  \n",
    "> ***–> To maximize the sample size for each regression, dropping NaN is performed after splitting.***  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "featureL = numDF.drop(columns=covarL).columns.tolist()\n",
    "\n",
    "#Split DF for each regression and drop NaN\n",
    "DF_splitL = []\n",
    "for feature in featureL:\n",
    "    tempL = [col_n for sublist in [[feature], covarL] for col_n in sublist]\n",
    "    tempDF = numDF[tempL].dropna()\n",
    "    DF_splitL.append(tempDF)\n",
    "\n",
    "print('The number of features:', len(DF_splitL))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Check examples for confirmation\n",
    "print(featureL[0])\n",
    "display(DF_splitL[0])\n",
    "print(featureL[1])\n",
    "display(DF_splitL[1])\n",
    "print(featureL[-1])\n",
    "display(DF_splitL[-1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Check each sample size\n",
    "tempL = []\n",
    "for feature_i in range(len(featureL)):\n",
    "    tempL.append(len(DF_splitL[feature_i]))\n",
    "tempS = pd.Series(tempL, index=featureL)\n",
    "tempS = tempS.sort_values(ascending=False)\n",
    "\n",
    "display(tempS)\n",
    "display(tempS.describe())\n",
    "\n",
    "#Distribution\n",
    "sns.set(style='ticks', font='Arial', context='notebook')\n",
    "plt.figure(figsize=(4, 3))\n",
    "sns.distplot(tempS)\n",
    "sns.despine()\n",
    "plt.xlabel('Sample size')\n",
    "plt.ylabel('Density')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1-3. Eliminate features with small sample size"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> No need for further elimination in this study cohort!  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1-4. Standardization of continuous numeric features and covariates"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Because all features are continuous numeric features this time, age is also standardized."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempL = []\n",
    "for feature_i in range(len(featureL)):\n",
    "    tempDF = DF_splitL[feature_i]\n",
    "    #Z-score transformation\n",
    "    tempDF1 = tempDF.select_dtypes(include=[np.number])\n",
    "    scaler = StandardScaler(copy=True, with_mean=True, with_std=True)\n",
    "    tempA = scaler.fit_transform(tempDF1)#Column direction\n",
    "    tempDF1 = pd.DataFrame(data=tempA, index=tempDF1.index, columns=tempDF1.columns)\n",
    "    ##Recover categorical covariates\n",
    "    tempDF2 = tempDF.select_dtypes(exclude=[np.number])\n",
    "    tempDF = pd.merge(tempDF1, tempDF2, left_index=True, right_index=True, how='left')\n",
    "    tempL.append(tempDF)\n",
    "DF_splitL = tempL#Update/overwrite"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Confirmation\n",
    "tempL = []\n",
    "for feature_i in range(len(featureL)):\n",
    "    tempL.append(DF_splitL[feature_i].loc[:, featureL[feature_i]])\n",
    "display(pd.concat(tempL, axis=1).describe())#Length is different but enforce merging to see the summary\n",
    "\n",
    "#Check distribution of some example numeric features\n",
    "sns.set(style='ticks', font='Arial', context='notebook')\n",
    "plt.figure(figsize=(4, 3))\n",
    "for feature_i in range(0, len(tempL), round(len(tempL)/4)):\n",
    "    sns.distplot(tempL[feature_i], label=tempL[feature_i].name)\n",
    "sns.despine()\n",
    "plt.xlabel('Z-score')\n",
    "plt.ylabel('Density')\n",
    "plt.legend(bbox_to_anchor=(1, 0.5), loc='center left', borderaxespad=1)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1-5. One-hot encoding for categorical covariates"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> category_encoders is more useful way than using sklearn.preprocessing or pandas.get_dummies.  \n",
    "> ***–> In statsmodels, categorical variables are automatically recognized!***  \n",
    "> –> Hence, this step is not needed anymore."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Check final DFs of examples for confirmation\n",
    "print(featureL[0])\n",
    "display(DF_splitL[0])\n",
    "print(featureL[1])\n",
    "display(DF_splitL[1])\n",
    "print(featureL[-1])\n",
    "display(DF_splitL[-1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Regression analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2-1. Perform OLS linear regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Model: (b)BMI ~ b0 + b1\\*Feature + b2\\*Sex + b3\\*Age + b4\\*AncestryPCs  \n",
    "> Main aim: Independently find associated features with BMI or biological BMI while adjusting sex, age and ancestry PCs as covariates  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempL = ['BMI', 'MetBMI', 'ProtBMI', 'ChemBMI', 'CombiBMI']\n",
    "\n",
    "#Initialize the result summary DF\n",
    "olsDF = pd.DataFrame(index=pd.Index(featureL, name='FeatureID'))\n",
    "olsDF = pd.merge(olsDF, numDF_meta['FeatureLabel'],\n",
    "                 left_index=True, right_index=True, how='left')\n",
    "olsDF['N'] = 1#Dummy int to keep the 2nd column, because it's same b/w BMI types in this cohort\n",
    "olsDF['DoF'] = 1#Dummy int to keep the 3rd column, because it's same b/w BMI types in this cohort\n",
    "\n",
    "#Perform OLS linear regression\n",
    "for bmi in tempL:\n",
    "    tempL1 = []#For R2 in the univariate model\n",
    "    tempL2 = []#For R2\n",
    "    tempL3 = []#For beta-coefficient\n",
    "    tempL4 = []#For SE of beta-coefficient\n",
    "    tempL5 = []#For 95% CI lower point of beta-coefficient (to use during visualization)\n",
    "    tempL6 = []#For 95% CI higher point of beta-coefficient (to use during visualization)\n",
    "    tempL7 = []#For t-statistic\n",
    "    tempL8 = []#For P-value\n",
    "    tempL9 = []#For sample size\n",
    "    tempL10 = []#For residual degrees of freedom\n",
    "    t_start = time.time()\n",
    "    for feature_i in range(len(featureL)):\n",
    "        feature = featureL[feature_i]\n",
    "        #Prepare DF\n",
    "        tempDF = DF_splitL[feature_i]\n",
    "        tempS = bmiDF['log_'+bmi]#Unstandardized\n",
    "        tempDF = pd.merge(tempDF, tempS, left_index=True, right_index=True, how='left')\n",
    "        #Rename dependent/independent variables\n",
    "        tempDF = tempDF.rename(columns={'log_'+bmi:'log_BMI', feature:'Feature'})\n",
    "        ##Add a constant for the intercept -> Similar to R, smf automatically add a constant\n",
    "        \n",
    "        #Fit univariate model\n",
    "        formula = 'log_BMI ~ Feature'\n",
    "        fit_res = smf.ols(formula, data=tempDF).fit()\n",
    "        #Save R2 [%]\n",
    "        tempL1.append(fit_res.rsquared*100)\n",
    "        \n",
    "        #Fit full model\n",
    "        formula = 'log_BMI ~ Feature + C(Sex) + Age + PC1 + PC2 + PC3 + PC4 + PC5'\n",
    "        fit_res = smf.ols(formula, data=tempDF).fit()\n",
    "        #Save R2 [%]\n",
    "        tempL2.append(fit_res.rsquared*100)\n",
    "        #Save beta-coefficient of the variable\n",
    "        tempL3.append(fit_res.params['Feature'])\n",
    "        tempL4.append(fit_res.bse['Feature'])\n",
    "        tempL5.append(fit_res.conf_int(alpha=0.05).loc['Feature', 0])\n",
    "        tempL6.append(fit_res.conf_int(alpha=0.05).loc['Feature', 1])\n",
    "        #Save t-statistic of the variable\n",
    "        tempL7.append(fit_res.tvalues['Feature'])\n",
    "        #Save P-value of the variable\n",
    "        tempL8.append(fit_res.pvalues['Feature'])\n",
    "        \n",
    "        #Save the common values b/w BMI types\n",
    "        tempL9.append(len(tempDF))#Same b/w BMI types in this cohort\n",
    "        tempL10.append(int(fit_res.df_resid))#Same b/w BMI types in this cohort\n",
    "    t_elapsed = time.time() - t_start\n",
    "    print('Elapsed time for', len(featureL), 'OLS linear regressions of log_'+bmi+':',\n",
    "          round(t_elapsed//60), 'min', round(t_elapsed%60, 1), 'sec')\n",
    "    \n",
    "    #Clean the results\n",
    "    tempDF = pd.DataFrame({bmi+'_UnivarR2':tempL1, bmi+'_R2':tempL2,\n",
    "                           bmi+'_Bcoef':tempL3, bmi+'_BcoefSE':tempL4,\n",
    "                           bmi+'_BcoefCIlow':tempL5, bmi+'_BcoefCIhigh':tempL6,\n",
    "                           bmi+'_tStat':tempL7, bmi+'_Pval':tempL8},\n",
    "                          index=pd.Index(featureL, name='FeatureID'))\n",
    "    ##P-value adjustment (within each BMI type) by using Benjamini–Hochberg method\n",
    "    tempDF[bmi+'_AdjPval_within'] = multi.multipletests(tempDF[bmi+'_Pval'], alpha=0.05, method='fdr_bh',\n",
    "                                                        is_sorted=False, returnsorted=False)[1]\n",
    "    ##Add dummy column for now\n",
    "    tempDF[bmi+'_AdjPval_all'] = 1.0#Dummy value\n",
    "    \n",
    "    #Merge the cleaned results\n",
    "    olsDF = pd.merge(olsDF, tempDF, left_index=True, right_index=True)\n",
    "\n",
    "#P-value adjustment (across all tests) by using Benjamini–Hochberg method\n",
    "tempL1 = olsDF.loc[:, olsDF.columns.str.contains('_Pval')].columns.tolist()\n",
    "tempDF = olsDF.reset_index().melt(var_name='BMItype', value_name='Pval', value_vars=tempL1,\n",
    "                                  id_vars=['FeatureID'])\n",
    "tempDF['AdjPval'] = multi.multipletests(tempDF['Pval'], alpha=0.05, method='fdr_bh',\n",
    "                                        is_sorted=False, returnsorted=False)[1]\n",
    "tempDF = tempDF.pivot(index='FeatureID', columns='BMItype', values='AdjPval')\n",
    "tempDF.columns =tempDF.columns.str.replace('_Pval', '_AdjPval_all')\n",
    "##Replace the dummy values with the adjusted p-values\n",
    "tempL1 = olsDF.loc[:, olsDF.columns.str.contains('_AdjPval_all')].columns.tolist()\n",
    "for col_n in tempL1:\n",
    "    olsDF[col_n] = tempDF[col_n]\n",
    "\n",
    "#Clean the result summary table\n",
    "olsDF['N'] = tempL9#Replace the dummy int values with the last values during the for-loop\n",
    "olsDF['DoF'] = tempL10#Replace the dummy int values with the last values during the for-loop\n",
    "olsDF = olsDF.sort_values(by=['BMI_Pval'], ascending=True)\n",
    "\n",
    "display(olsDF)\n",
    "#Save\n",
    "fileDir = './ExportData/'\n",
    "ipynbName = '220803_Multiomics-BMI-NatMed1stRevision_BMI-regression_NumericFeatures_'\n",
    "fileName = 'result-summary.tsv'\n",
    "olsDF.to_csv(fileDir+ipynbName+fileName, sep='\\t', index=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2-2. Significantly associated features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Significantly associated features\n",
    "print('Significantly associated features (FDR < 0.05):')\n",
    "tempD = {'BMI':'0.3', 'MetBMI':'b', 'ProtBMI':'r', 'ChemBMI':'g', 'CombiBMI':'m'}\n",
    "tempL = []\n",
    "for bmi in list(tempD.keys()):\n",
    "    #Extact significantly associated features\n",
    "    tempDF = olsDF.loc[olsDF[bmi+'_AdjPval_all']<0.05]\n",
    "    tempL.append(tempDF.index.tolist())\n",
    "    \n",
    "    #Cf. Adjustment within each BMI type\n",
    "    tempDF1 = olsDF.loc[olsDF[bmi+'_AdjPval_within']<0.05]\n",
    "    print(' - '+bmi+':', len(tempDF), '(cf.', len(tempDF1), 'by the within-adjustment)')\n",
    "##Flatten and drop multiplicates\n",
    "tempL = list(set(item for sublist in tempL for item in sublist))\n",
    "tempDF = olsDF.loc[olsDF.index.isin(tempL)]\n",
    "display(tempDF)\n",
    "\n",
    "#Add sample size to feature label\n",
    "tempDF = tempDF.reset_index()\n",
    "tempDF['N'] = [f'{item:,}' for item in tempDF['N'].tolist()]\n",
    "tempDF['FeatureID'] = tempDF['FeatureID'].str.cat(tempDF['N'], sep=' ('+r'$n$'+' = ')\n",
    "tempDF['FeatureID'] = tempDF['FeatureID'].str.cat(np.repeat(')', len(tempDF)), sep='')\n",
    "tempDF = tempDF.set_index('FeatureID')\n",
    "\n",
    "#Visualize R2\n",
    "tempDF1 = tempDF[[item+'_R2' for item in list(tempD.keys())]].sort_values(by=['BMI_R2'], ascending=False)\n",
    "tempL = tempDF1.index.tolist()#Save order\n",
    "tempDF1 = tempDF1.reset_index().melt(var_name='BMItype', value_name='R2', id_vars=['FeatureID'])\n",
    "tempDF1['BMItype'] = tempDF1['BMItype'].str.replace('_R2', '')\n",
    "##Style and annotation info\n",
    "tempL1 = []\n",
    "tempL2 = []\n",
    "for row_i in range(len(tempDF1)):\n",
    "    bmi = tempDF1.iloc[row_i]['BMItype']\n",
    "    feature = tempDF1.iloc[row_i]['FeatureID']\n",
    "    #Association\n",
    "    if tempDF.loc[feature, bmi+'_Bcoef']>0:\n",
    "        tempL1.append('Positive association')\n",
    "    elif tempDF.loc[feature, bmi+'_Bcoef']<0:\n",
    "        tempL1.append('Negative association')\n",
    "    else:#just in case\n",
    "        tempL1.append('No association')\n",
    "    #Significance\n",
    "    if tempDF.loc[feature, bmi+'_AdjPval_all']<0.05:\n",
    "        tempL2.append('FDR<0.05')\n",
    "    else:\n",
    "        tempL2.append('n.s.')\n",
    "tempDF1['Association'] = tempL1\n",
    "tempDF1['Signif'] = tempL2\n",
    "##Plot\n",
    "sns.set(style='ticks', font='Arial', context='talk')\n",
    "plt.figure(figsize=(7.5, 18))\n",
    "p = sns.barplot(data=tempDF1, x='R2', y='FeatureID', order=tempL, hue='BMItype', hue_order=tempD.keys(),\n",
    "                palette=tempD, dodge=True, edgecolor='black', linewidth=1)\n",
    "p.grid(axis='x', linestyle='--', color='black')\n",
    "for rect, row_i in zip(p.patches, range(len(tempDF1))):\n",
    "    xcoord = rect.get_width()\n",
    "    ycoord = rect.get_y() + rect.get_height()/2\n",
    "    #Add annotation\n",
    "    label = tempDF1.iloc[row_i]['Signif']\n",
    "    bmi = tempDF1.iloc[row_i]['BMItype']\n",
    "    if label=='n.s.':\n",
    "        text_offset = 50\n",
    "        hue_order = list(tempD.keys()).index(bmi)\n",
    "        if xcoord > 0:\n",
    "            if hue_order%2 == 0:\n",
    "                offset = +text_offset/2\n",
    "            else:\n",
    "                offset = +text_offset\n",
    "            halign = 'left'\n",
    "        else:\n",
    "            if hue_order%2 == 0:\n",
    "                offset = -text_offset/2\n",
    "            else:\n",
    "                offset = -text_offset\n",
    "            halign = 'right'\n",
    "        p.annotate(label, xy=(xcoord, ycoord), xytext=(offset, 0), textcoords='offset points',\n",
    "                   horizontalalignment=halign, verticalalignment='center',\n",
    "                   fontsize='x-small', color='black',\n",
    "                   arrowprops={'arrowstyle':'-', 'color':tempD[bmi], 'linewidth':1,\n",
    "                               'shrinkA':offset/25, 'shrinkB':offset/5})\n",
    "sns.despine()\n",
    "plt.xlabel('Ratio of explained variance [%]')\n",
    "plt.ylabel('')\n",
    "for ycoord in range(len(tempDF)):\n",
    "    if ycoord%2 == 0:\n",
    "        plt.axhspan(ymin=ycoord-0.5, ymax=ycoord+0.5, facecolor='k', alpha=0.2)\n",
    "plt.margins(0.02, 0.005, tight=True)\n",
    "plt.legend(loc='lower right')\n",
    "plt.show()\n",
    "\n",
    "#Visualize beta-coefficient\n",
    "tempDF1 = tempDF[[item+'_Bcoef' for item in list(tempD.keys())]].sort_values(by=['BMI_Bcoef'], ascending=False)\n",
    "tempL = tempDF1.index.tolist()#Save order\n",
    "tempDF1 = tempDF1.reset_index().melt(var_name='BMItype', value_name='Bcoef', id_vars=['FeatureID'])\n",
    "tempDF1['BMItype'] = tempDF1['BMItype'].str.replace('_Bcoef', '')\n",
    "##Annotation info\n",
    "tempL1 = []\n",
    "tempL2 = []\n",
    "for row_i in range(len(tempDF1)):\n",
    "    bmi = tempDF1.iloc[row_i]['BMItype']\n",
    "    feature = tempDF1.iloc[row_i]['FeatureID']\n",
    "    #SE\n",
    "    tempL1.append(tempDF.loc[feature, bmi+'_BcoefSE'])\n",
    "    #Significance\n",
    "    if tempDF.loc[feature, bmi+'_AdjPval_all']<0.05:\n",
    "        tempL2.append('FDR<0.05')\n",
    "    else:\n",
    "        tempL2.append('n.s.')\n",
    "tempDF1['BcoefSE'] = tempL1\n",
    "tempDF1['Signif'] = tempL2\n",
    "##Plot\n",
    "sns.set(style='ticks', font='Arial', context='talk')\n",
    "plt.figure(figsize=(7.5, 18))\n",
    "p = sns.barplot(data=tempDF1, x='Bcoef', y='FeatureID', order=tempL, hue='BMItype', hue_order=tempD.keys(),\n",
    "                palette=tempD, dodge=True, edgecolor='white', linewidth=0.5, saturation=1)\n",
    "p.grid(axis='x', linestyle='--', color='black')\n",
    "for rect, row_i in zip(p.patches, range(len(tempDF1))):\n",
    "    xcoord = rect.get_width()\n",
    "    ycoord = rect.get_y() + rect.get_height()/2\n",
    "    #Add SEM\n",
    "    sem = tempDF1.iloc[row_i]['BcoefSE']\n",
    "    plt.errorbar(x=xcoord, y=ycoord, xerr=sem, fmt='none', ecolor='k', elinewidth=1, capsize=2)\n",
    "    #Add annotation\n",
    "    label = tempDF1.iloc[row_i]['Signif']\n",
    "    bmi = tempDF1.iloc[row_i]['BMItype']\n",
    "    if label=='n.s.':\n",
    "        text_offset = 50\n",
    "        hue_order = list(tempD.keys()).index(bmi)\n",
    "        if xcoord > 0:\n",
    "            if hue_order%2 == 0:\n",
    "                offset = +text_offset/2\n",
    "            else:\n",
    "                offset = +text_offset\n",
    "            halign = 'left'\n",
    "            xcoord = xcoord + sem\n",
    "        else:\n",
    "            if hue_order%2 == 0:\n",
    "                offset = -text_offset/2\n",
    "            else:\n",
    "                offset = -text_offset\n",
    "            halign = 'right'\n",
    "            xcoord = xcoord - sem\n",
    "        p.annotate(label, xy=(xcoord, ycoord), xytext=(offset, 0), textcoords='offset points',\n",
    "                   horizontalalignment=halign, verticalalignment='center',\n",
    "                   fontsize='x-small', color='black',\n",
    "                   arrowprops={'arrowstyle':'-', 'color':tempD[bmi], 'linewidth':1,\n",
    "                               'shrinkA':offset/25, 'shrinkB':offset/5})\n",
    "sns.despine()\n",
    "plt.xlabel(r'$\\beta$'+'-coefficient')\n",
    "plt.ylabel('')\n",
    "for ycoord in range(len(tempDF)):\n",
    "    if ycoord%2 == 0:\n",
    "        plt.axhspan(ymin=ycoord-0.5, ymax=ycoord+0.5, facecolor='k', alpha=0.2)\n",
    "plt.margins(0.02, 0.005, tight=True)\n",
    "plt.legend(loc='lower right')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2-3. For paper figure"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Basically, P < 0.05 should be described. However, because of tiny figure, it is abbreviated and the exact P-value will be provided as a supplementary table.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Significantly associated features\n",
    "print('Significantly associated features (FDR < 0.05):')\n",
    "tempD = {'BMI':'0.3', 'MetBMI':'b', 'ProtBMI':'r', 'ChemBMI':'g', 'CombiBMI':'m'}\n",
    "tempL = []\n",
    "for bmi in list(tempD.keys()):\n",
    "    #Extact significantly associated features\n",
    "    tempDF = olsDF.loc[olsDF[bmi+'_AdjPval_all']<0.05]\n",
    "    tempL.append(tempDF.index.tolist())\n",
    "    print(' - '+bmi+':', len(tempDF))\n",
    "##Flatten and drop multiplicates\n",
    "tempL = list(set(item for sublist in tempL for item in sublist))\n",
    "tempDF = olsDF.loc[olsDF.index.isin(tempL)]\n",
    "\n",
    "#Add sample size to feature label\n",
    "tempDF = tempDF.reset_index()\n",
    "tempDF['N'] = [f'{item:,}' for item in tempDF['N'].tolist()]\n",
    "tempDF['FeatureLabel'] = tempDF['FeatureLabel'].str.cat(tempDF['N'], sep='\\n('+r'$n$'+' = ')\n",
    "tempDF['FeatureLabel'] = tempDF['FeatureLabel'].str.cat(np.repeat(')', len(tempDF)), sep='')\n",
    "tempDF = tempDF.set_index('FeatureLabel')\n",
    "\n",
    "#Visualize beta-coefficient with forest plot\n",
    "##Sort feature\n",
    "tempDF1 = tempDF[[item+'_Bcoef' for item in tempD.keys()]].sort_values(by=['BMI_Bcoef'], ascending=False)\n",
    "tempL = tempDF1.index.tolist()#Save order\n",
    "##Clean DF for plot\n",
    "tempDF1 = tempDF.reset_index().melt(var_name='BMItype', value_name='Bcoef', id_vars=['FeatureLabel'],\n",
    "                                    value_vars=[item+'_Bcoef' for item in tempD.keys()])\n",
    "tempDF1['BMItype'] = tempDF1['BMItype'].str.replace('_Bcoef', '')\n",
    "tempDF2 = tempDF.reset_index().melt(var_name='BMItype', value_name='BcoefCIlow', id_vars=['FeatureLabel'],\n",
    "                                    value_vars=[item+'_BcoefCIlow' for item in tempD.keys()])\n",
    "tempDF2['BMItype'] = tempDF2['BMItype'].str.replace('_BcoefCIlow', '')\n",
    "tempDF1 = pd.merge(tempDF1, tempDF2, on=['FeatureLabel', 'BMItype'], how='left')\n",
    "tempDF2 = tempDF.reset_index().melt(var_name='BMItype', value_name='BcoefCIhigh', id_vars=['FeatureLabel'],\n",
    "                                    value_vars=[item+'_BcoefCIhigh' for item in tempD.keys()])\n",
    "tempDF2['BMItype'] = tempDF2['BMItype'].str.replace('_BcoefCIhigh', '')\n",
    "tempDF1 = pd.merge(tempDF1, tempDF2, on=['FeatureLabel', 'BMItype'], how='left')\n",
    "##Convert CI points to difference values\n",
    "tempDF1['BcoefCIlow'] = tempDF1['Bcoef'] - tempDF1['BcoefCIlow']\n",
    "tempDF1['BcoefCIhigh'] = tempDF1['BcoefCIhigh'] - tempDF1['Bcoef']\n",
    "##Significance annotation\n",
    "tempL1 = []\n",
    "for row_i in range(len(tempDF1)):\n",
    "    bmi = tempDF1['BMItype'].iloc[row_i]\n",
    "    feature = tempDF1['FeatureLabel'].iloc[row_i]\n",
    "    #Significance\n",
    "    if tempDF.loc[feature, bmi+'_AdjPval_all']<0.001:\n",
    "        tempL1.append('***')\n",
    "    elif tempDF.loc[feature, bmi+'_AdjPval_all']<0.01:\n",
    "        tempL1.append('**')\n",
    "    elif tempDF.loc[feature, bmi+'_AdjPval_all']<0.05:\n",
    "        tempL1.append('*')\n",
    "    else:\n",
    "        tempL1.append('n.s.')\n",
    "tempDF1['Signif'] = tempL1\n",
    "##Color for edge\n",
    "tempDF1['EdgeColor'] = tempDF1['BMItype'].map(tempD)\n",
    "##Plot\n",
    "ncols = 2\n",
    "sns.set(style='ticks', font='Arial', context='talk')\n",
    "fig, axes = plt.subplots(nrows=1, ncols=ncols, figsize=(14, 11), sharex=False, sharey=False)\n",
    "for ax_i, ax in enumerate(axes.flat):\n",
    "    half = int(len(tempL)/ncols)\n",
    "    tempL1 = tempL[half*ax_i:half*(ax_i+1)]\n",
    "    tempDF2 = tempDF1.loc[tempDF1['FeatureLabel'].isin(tempL1)].reset_index()\n",
    "    #Point\n",
    "    sns.pointplot(data=tempDF2, x='Bcoef', y='FeatureLabel', order=tempL1,\n",
    "                  hue='BMItype', hue_order=tempD.keys(), palette=tempD,\n",
    "                  dodge=0.7, join=False, ci=None, markers='o', scale=0.6, ax=ax)\n",
    "    #Get coordinate of each point\n",
    "    xcoordL = []\n",
    "    ycoordL = []\n",
    "    for coord in ax.collections:\n",
    "        for x, y in coord.get_offsets():\n",
    "            xcoordL.append(x)\n",
    "            ycoordL.append(y)\n",
    "    #Add errorbars and annotation\n",
    "    for point_i in range(len(tempDF2)):\n",
    "        #Be careful about order of ax.collections: tempL -> tempD.keys()\n",
    "        label = tempL1[point_i%len(tempL1)]\n",
    "        bmi = list(tempD.keys())[point_i//len(tempL1)]\n",
    "        row_i = tempDF2.loc[(tempDF2['FeatureLabel']==label)&(tempDF2['BMItype']==bmi)].index.tolist()[0]\n",
    "        #Add errorbars manually\n",
    "        ax.errorbar(x=xcoordL[point_i], y=ycoordL[point_i],\n",
    "                    xerr=[[tempDF2['BcoefCIlow'].iloc[row_i]], [tempDF2['BcoefCIhigh'].iloc[row_i]]],\n",
    "                    fmt='', ecolor=tempDF2['EdgeColor'].iloc[row_i], elinewidth=4, capsize=0, capthick=4,\n",
    "                    linestyle='', zorder=1)\n",
    "        #Add significance\n",
    "        signif = tempDF2['Signif'].iloc[row_i]\n",
    "        if signif!='n.s.':\n",
    "            text_offset = 5\n",
    "            if xcoordL[point_i] > 0:\n",
    "                xoffset = +text_offset\n",
    "                halign = 'left'\n",
    "                xcoord = xcoordL[point_i] + tempDF2['BcoefCIhigh'].iloc[row_i]\n",
    "            else:\n",
    "                xoffset = -text_offset\n",
    "                halign = 'right'\n",
    "                xcoord = xcoordL[point_i] - tempDF2['BcoefCIlow'].iloc[row_i]\n",
    "            ycoord = ycoordL[point_i]\n",
    "            yoffset = -4#Because asterisk looks upper shift\n",
    "            ax.annotate(signif, xy=(xcoord, ycoord), xytext=(xoffset, yoffset), textcoords='offset points',\n",
    "                        horizontalalignment=halign, verticalalignment='center',\n",
    "                        fontsize='medium', color=tempD[bmi], zorder=2)\n",
    "    #Add reference line\n",
    "    #p.grid(axis='x', linestyle='--', color='k')\n",
    "    ax.axvline(x=0, **{'linestyle':'--', 'color':'k', 'zorder':0})\n",
    "    #Add shading\n",
    "    for ycoord in range(len(tempL1)):\n",
    "        if ycoord%2 == 0:\n",
    "            ax.axhspan(ymin=ycoord-0.5, ymax=ycoord+0.5, facecolor='k', alpha=0.2, zorder=0)\n",
    "    #Ax-dependent setting\n",
    "    if ax_i == 0:\n",
    "        #Range\n",
    "        ax.set(xlim=(-0.02, 0.22), xticks=np.arange(0.0, 0.21, 0.1))\n",
    "        #Add legend\n",
    "        ax.get_legend().remove()\n",
    "    elif ax_i ==1:\n",
    "        #Range\n",
    "        ax.set(xlim=(-0.12, 0.12), xticks=np.arange(-0.1, 0.11, 0.1))\n",
    "        #Add legend\n",
    "        ax.legend(title='BMI type', loc='lower right', handlelength=1.0, handletextpad=0.5)\n",
    "sns.despine()\n",
    "plt.setp(axes, ylim=(len(tempL1)-0.5, -0.5))#Otherwise, axes are extended\n",
    "plt.setp(axes, xlabel='', ylabel='')\n",
    "fig.tight_layout()\n",
    "fig.text(x=0.6175, y=0.0175,#Manual adjustment\n",
    "         s=r'$\\beta$'+'-coefficient [log-scaled measured or omics-inferred BMI (kg m'+r'$^{-2}$'+') per s.d.]',\n",
    "         fontsize='medium', verticalalignment='top', horizontalalignment='center')\n",
    "##Save\n",
    "fileDir = './ExportFigures/'\n",
    "ipynbName = '220803_Multiomics-BMI-NatMed1stRevision_BMI-regression_NumericFeatures_'\n",
    "fileName = 'bcoef.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": [
    "# — 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
}