Code23_WHtR-LASSO-bcoef.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Multiomics BMI Paper — Beta-coefficients of WHtR LASSO Models"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***by Kengo Watanabe***  \n",
    "\n",
    "This Jupyter Notebook (with Python 3 kernel) visualized beta-coefficients of the blood omics-based WHtR LASSO models, and assessed relationships between WHtR and the retained variables using regression analysis (in the Arivale cohort).  \n",
    "\n",
    "Input files:  \n",
    "* Arivale baseline WHtR: 220621_Multiomics-BMI-NatMedRevision_WHtR-DataCleaning_baseline-WHtR-final-cohort.tsv  \n",
    "* Arivale baseline blood omics (preprocessed): 210104_Biological-BMI-paper_RF-imputation_baseline-\\[metDF/protDF/chemDF/combiDF\\]-with-RF-imputation.tsv  \n",
    "* Arivale baseline WHtR predictions: 220822_Multiomics-BMI-NatMed1stRevision_WHtR-baseline-LASSO-ver2_\\[MetWHtR/ProtWHtR/ChemWHtR/CombiWHtR\\]-\\[Female/Male/BothSex\\].tsv  \n",
    "* Biological WHtR models: 220822_Multiomics-BMI-NatMed1stRevision_WHtR-baseline-LASSO-ver2_\\[MetWHtR/ProtWHtR/ChemWHtR/CombiWHtR\\]-\\[Female/Male/BothSex\\]-LASSObcoefs.tsv  \n",
    "\n",
    "Output figures and tables:  \n",
    "* Supplementary Figure 7j, 7k  \n",
    "* Tables for Supplementary Data 9  \n",
    "\n",
    "Original notebook (memo for my future tracing):  \n",
    "* dalek:\\[JupyterLab HOME\\]/220621_Multiomics-BMI-NatMedRevision/220823_Multiomics-BMI-NatMed1stRevision_WHtR-LASSO-bcoef-ver2.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",
    "#!conda install -c conda-forge -y matplotlib-venn\n",
    "from matplotlib_venn import venn3, venn3_circles, venn2, venn2_circles\n",
    "import statsmodels.formula.api as smf\n",
    "from statsmodels.stats import multitest as multi\n",
    "import matplotlib.patches as mpatches\n",
    "\n",
    "!conda list"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Data preparation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1-1. Import the cleaned dataframes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Import the baseline WHtR dataframe\n",
    "fileDir = './ExportData/'\n",
    "ipynbName = '220621_Multiomics-BMI-NatMedRevision_WHtR-DataCleaning_'\n",
    "fileName = 'baseline-WHtR-final-cohort.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",
    "##Take WHtR and general covariates (without Race in this study)\n",
    "tempL = ['log_BaseWHtR', 'log_BaseBMI', 'Sex', 'BaseAge', 'PC1', 'PC2', 'PC3', 'PC4', 'PC5']\n",
    "tempDF = tempDF[tempL]\n",
    "\n",
    "display(tempDF)\n",
    "\n",
    "whtrDF = tempDF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempDF1 = whtrDF\n",
    "\n",
    "#Import the cleaned baseline omics dataframes\n",
    "fileDir = '../210104_Biological-BMI-paper/ExportData/'\n",
    "ipynbName = '210104_Biological-BMI-paper_RF-imputation_'\n",
    "tempL = ['log_BaseBMI', 'Sex', 'BaseAge', 'PC1', 'PC2', 'PC3', 'PC4', 'PC5', 'Race']\n",
    "tempD = {}\n",
    "for df_n in ['metDF', 'chemDF', 'protDF', 'combiDF']:\n",
    "    fileName = 'baseline-'+df_n+'-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",
    "    print(df_n+' original shape:', tempDF.shape)\n",
    "    #Drop BMI and covariates\n",
    "    tempDF = tempDF.drop(columns=tempL)\n",
    "    #Extract the individuals having WHtR\n",
    "    tempDF = tempDF.loc[tempDF1.index.tolist()]\n",
    "    display(tempDF)\n",
    "    tempD[df_n] = tempDF\n",
    "\n",
    "metDF = tempD['metDF']\n",
    "chemDF = tempD['chemDF']\n",
    "protDF = tempD['protDF']\n",
    "combiDF = tempD['combiDF']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1-2. Stratification with sex"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Stratify the cohort with sex\n",
    "whtrDF_F = whtrDF.loc[whtrDF['Sex']=='F']\n",
    "whtrDF_M = whtrDF.loc[whtrDF['Sex']=='M']\n",
    "whtrDF_B = whtrDF#Not copy just rename\n",
    "print('Female, Male, Both sex = ', len(whtrDF_F), ', ', len(whtrDF_M), ', ', len(whtrDF_B))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Metabolomics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2-1. Prepare and standardize variables"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Variables, including bWHtR (as a reference) and covariates, are standardized for OLS linear regression.  \n",
    "> –> In this version, unstandardized WHtR is used as the dependent variable. However, WHtR within the following DFs is standardized (just for the code simplicity).  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempDF = metDF\n",
    "tempD1 = {'Female':whtrDF_F, 'Male':whtrDF_M, 'BothSex':whtrDF_B}\n",
    "bwhtr = 'MetWHtR'\n",
    "tempD2 = {}\n",
    "for sex in tempD1.keys():\n",
    "    #Add covariates (and WHtR) while selecting the cohort\n",
    "    tempDF1 = tempD1[sex]\n",
    "    tempDF1 = pd.merge(tempDF, tempDF1, left_index=True, right_index=True, how='inner')\n",
    "    \n",
    "    #Add biological WHtR\n",
    "    fileDir = './ExportData/'\n",
    "    ipynbName = '220822_Multiomics-BMI-NatMed1stRevision_WHtR-baseline-LASSO-ver2_'\n",
    "    fileName = bwhtr+'-'+sex+'.tsv'\n",
    "    tempDF2 = pd.read_csv(fileDir+ipynbName+fileName, sep='\\t', dtype={'public_client_id':str})\n",
    "    tempDF2 = tempDF2.set_index('public_client_id')\n",
    "    tempDF1['log_Base'+bwhtr] = tempDF2['log_Base'+bwhtr]\n",
    "    \n",
    "    #Z-score transformation\n",
    "    tempDF2 = tempDF1.select_dtypes(include=[np.number])\n",
    "    scaler = StandardScaler(copy=True, with_mean=True, with_std=True)\n",
    "    tempA = scaler.fit_transform(tempDF2)#Column direction\n",
    "    tempDF2 = pd.DataFrame(data=tempA, index=tempDF2.index, columns=tempDF2.columns)\n",
    "    ##Recover categorical covariates\n",
    "    tempDF3 = tempDF1.select_dtypes(exclude=[np.number])\n",
    "    tempDF1 = pd.merge(tempDF2, tempDF3, left_index=True, right_index=True, how='left')\n",
    "    \n",
    "    print(sex)\n",
    "    display(tempDF1.describe(include='all'))\n",
    "    tempD2[sex] = tempDF1\n",
    "\n",
    "metDF_F = tempD2['Female']\n",
    "metDF_M = tempD2['Male']\n",
    "metDF_B = tempD2['BothSex']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2-2. Correlation b/w all pairwise variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Female and male cohorts\n",
    "tempD1 = {'Female':metDF_F, 'Male':metDF_M}\n",
    "bwhtr = 'MetWHtR'\n",
    "tempD2 = {}\n",
    "for sex in tempD1.keys():\n",
    "    print(sex+' cohort:')\n",
    "    tempDF = tempD1[sex]\n",
    "    #Compute correlation matrix\n",
    "    tempL = [item for sublist in [whtrDF.columns.tolist(), ['log_Base'+bwhtr]] for item in sublist]\n",
    "    tempDF1 = tempDF.drop(columns=tempL)\n",
    "    tempDF1 = tempDF1.corr(method='pearson')\n",
    "    print('• Combinations:', int(len(tempDF1)*(len(tempDF1)-1)/2))\n",
    "    \n",
    "    #Extract lower triangle matrix\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",
    "    print('• nrows:', len(tempDF1))\n",
    "    print('• |Pearson\\'s r| > 0.8:', len(tempDF1.loc[abs(tempDF1['Pearson_r'])>0.8]))\n",
    "    display(tempDF1.loc[abs(tempDF1['Pearson_r'])>0.8])\n",
    "    \n",
    "    tempD2[sex] = tempDF1\n",
    "\n",
    "#Distribution\n",
    "tempD3 = {'Female':'tab:red', 'Male':'tab:blue'}\n",
    "sns.set(style='ticks', font='Arial', context='notebook')\n",
    "plt.figure(figsize=(4, 3))\n",
    "for sex in tempD2.keys():\n",
    "    tempS = tempD2[sex]['Pearson_r']\n",
    "    sns.distplot(tempS, color=tempD3[sex], label=sex).set(xlim=(-1, 1), xticks=np.arange(-1, 1.1, 0.5))\n",
    "sns.despine()\n",
    "plt.xlabel('Pearson\\'s '+r'$r$')\n",
    "plt.ylabel('Density')\n",
    "plt.legend(loc='upper right')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Both sex cohort\n",
    "tempDF = metDF_B\n",
    "bwhtr = 'MetWHtR'\n",
    "bwhtr_color = 'b'\n",
    "\n",
    "print('Both sex cohort:')\n",
    "#Compute correlation matrix\n",
    "tempL = [item for sublist in [whtrDF.columns.tolist(), ['log_Base'+bwhtr]] for item in sublist]\n",
    "tempDF = tempDF.drop(columns=tempL)\n",
    "tempDF = tempDF.corr(method='pearson')\n",
    "print('• Combinations:', int(len(tempDF)*(len(tempDF)-1)/2))\n",
    "\n",
    "#Extract lower triangle matrix\n",
    "tempDF = tempDF.where(np.tril(np.ones(tempDF.shape), k=-1).astype(np.bool), other=np.nan)\n",
    "tempDF.index.rename('Variable1', inplace=True)\n",
    "tempDF = tempDF.reset_index().melt(var_name='Variable2', value_name='Pearson_r', id_vars=['Variable1'])\n",
    "tempDF = tempDF.dropna()\n",
    "print('• nrows:', len(tempDF))\n",
    "print('• |Pearson\\'s r| > 0.8:', len(tempDF.loc[abs(tempDF['Pearson_r'])>0.8]))\n",
    "display(tempDF.loc[abs(tempDF['Pearson_r'])>0.8])\n",
    "\n",
    "#Distribution\n",
    "sns.set(style='ticks', font='Arial', context='talk')\n",
    "plt.figure(figsize=(4, 3))\n",
    "sns.distplot(tempDF['Pearson_r'], color=bwhtr_color).set(xlim=(-1, 1), xticks=np.arange(-1, 1.1, 0.5))\n",
    "plt.axvline(x=tempDF['Pearson_r'].max(), **{'linestyle':'--', 'color':bwhtr_color})\n",
    "plt.axvline(x=tempDF['Pearson_r'].min(), **{'linestyle':'--', 'color':bwhtr_color})\n",
    "sns.despine()\n",
    "plt.ylabel('Density')\n",
    "plt.xlabel('Pearson\\'s '+r'$r$')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2-3. Prepare the beta-coefficients of LASSO models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bwhtr = 'MetWHtR'\n",
    "tempD = {}\n",
    "for sex in ['Female', 'Male', 'BothSex']:\n",
    "    #Import the LASSO beta-coefficients\n",
    "    fileDir = './ExportData/'\n",
    "    ipynbName = '220822_Multiomics-BMI-NatMed1stRevision_WHtR-baseline-LASSO-ver2_'\n",
    "    fileName = bwhtr+'-'+sex+'-LASSObcoefs.tsv'\n",
    "    tempDF = pd.read_csv(fileDir+ipynbName+fileName, sep='\\t').set_index('Variable')\n",
    "    tempDF = tempDF.drop(index=['Intercept'])\n",
    "    tempD[sex] = tempDF\n",
    "    \n",
    "    #Check\n",
    "    print(sex+':')\n",
    "    print(' - Variables:', len(tempDF))\n",
    "    tempDF1 = tempDF.loc[tempDF['nZeros']!=10]\n",
    "    print(' - Variables with non-zero beta-coefficient:', len(tempDF1),\n",
    "          '(', len(tempDF1)/len(tempDF)*100, '%)')\n",
    "    tempDF1 = tempDF.loc[tempDF['nZeros']==0]\n",
    "    print(' - Variables with non-zero beta-coefficient in all 10 models:', len(tempDF1),\n",
    "          '(', len(tempDF1)/len(tempDF)*100, '%)')\n",
    "    tempDF1 = tempDF1.sort_values(by='Mean', ascending=False)\n",
    "    display(tempDF1)\n",
    "    print('')\n",
    "\n",
    "metWHtR_F_bcoefs = tempD['Female']\n",
    "metWHtR_M_bcoefs = tempD['Male']\n",
    "metWHtR_B_bcoefs = tempD['BothSex']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2-4. Variables retained in all 10 models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempDF1 = metWHtR_F_bcoefs\n",
    "tempDF2 = metWHtR_M_bcoefs\n",
    "tempDF3 = metWHtR_B_bcoefs\n",
    "category = 'Metabolomics'\n",
    "\n",
    "#Variables with non-zero beta-coefficient in all 10 models\n",
    "tempS1 = set(tempDF1.loc[tempDF1['nZeros']==0].index)\n",
    "tempS2 = set(tempDF2.loc[tempDF2['nZeros']==0].index)\n",
    "tempS3 = set(tempDF3.loc[tempDF3['nZeros']==0].index)\n",
    "\n",
    "#Common the robust variables\n",
    "var_111 = tempS1 & tempS2 & tempS3\n",
    "print('Common variables with non-zero beta-coefficient in all 10 models:\\n', var_111)\n",
    "var_110 = (tempS1 & tempS2) - var_111\n",
    "var_101 = (tempS1 & tempS3) - var_111\n",
    "var_011 = (tempS2 & tempS3) - var_111\n",
    "var_100 = tempS1 - (var_110 | var_101 | var_111)\n",
    "var_010 = tempS2 - (var_110 | var_011 | var_111)\n",
    "var_001 = tempS3 - (var_101 | var_011 | var_111)\n",
    "\n",
    "#The number of each subset\n",
    "var_subset = [len(var_100), len(var_010), len(var_110), len(var_001),\n",
    "              len(var_101), len(var_011), len(var_111)]\n",
    "print('Each subset:', var_subset)\n",
    "\n",
    "#Venn diagram\n",
    "sns.set(font='Arial', context='talk')\n",
    "venn3(subsets=var_subset, set_labels=('Female', 'Male', 'Both sex'),\n",
    "      set_colors=('tab:red', 'tab:blue', 'tab:green'), alpha=0.4)\n",
    "venn3_circles(subsets=var_subset)\n",
    "plt.title('Robust variables\\n—'+category+'—', fontdict={'fontsize':24})\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempDF1 = metWHtR_F_bcoefs\n",
    "tempDF2 = metWHtR_M_bcoefs\n",
    "tempDF3 = metWHtR_B_bcoefs\n",
    "bwhtr_color = 'b'\n",
    "\n",
    "#Variables with non-zero beta-coefficient in all 10 models\n",
    "tempS1 = set(tempDF1.loc[tempDF1['nZeros']==0].index)\n",
    "tempS2 = set(tempDF2.loc[tempDF2['nZeros']==0].index)\n",
    "tempS3 = set(tempDF3.loc[tempDF3['nZeros']==0].index)\n",
    "\n",
    "#Union of the robust variables\n",
    "tempL = list(tempS1 | tempS2 | tempS3)\n",
    "print('Union of the variables with non-zero beta-coefficient in all 10 models:', len(tempL))\n",
    "\n",
    "#Merge (use original bcoefs to rescue ones with zeros in any 10 models)\n",
    "tempDF1 = tempDF1.loc[tempDF1.index.isin(tempL)].drop(columns=['Mean', 'SD', 'nZeros'])\n",
    "tempDF2 = tempDF2.loc[tempDF2.index.isin(tempL)].drop(columns=['Mean', 'SD', 'nZeros'])\n",
    "tempDF3 = tempDF3.loc[tempDF3.index.isin(tempL)].drop(columns=['Mean', 'SD', 'nZeros'])\n",
    "tempDF = pd.concat([tempDF1, tempDF2, tempDF3], axis=0)\n",
    "tempDF['Model'] = np.repeat(['Female', 'Male', 'Both sex'], [len(tempDF1), len(tempDF2), len(tempDF3)])\n",
    "\n",
    "#tidyr::gather in R\n",
    "tempDF1 = tempDF.reset_index().melt(var_name='LASSO', value_name='bcoef', id_vars=['Variable', 'Model'])\n",
    "\n",
    "#Order for plot\n",
    "tempDF2 = tempDF1.loc[tempDF1['Model']=='Both sex'].groupby('Variable', as_index=True).agg({'bcoef':np.mean})\n",
    "tempDF2 = tempDF2.sort_values(by='bcoef', ascending=False)\n",
    "\n",
    "#Plot\n",
    "sns.set(style='ticks', font='Arial', context='talk')\n",
    "sns.catplot(data=tempDF1, y='Variable', order=tempDF2.index.tolist(), x='bcoef', kind='box',\n",
    "            hue='Model', palette='Set1', dodge=True, height=35, aspect=0.35,\n",
    "            showfliers=True, flierprops={'marker':'o', 'markerfacecolor':'gray', 'alpha':0.4})\n",
    "plt.grid(axis='x', linestyle='--', color='black')\n",
    "for row_i in range(len(tempDF2)):\n",
    "    if row_i%2 == 0:\n",
    "        plt.axhspan(ymin=row_i-0.5, ymax=row_i+0.5, facecolor=bwhtr_color, alpha=0.2, zorder=0)\n",
    "plt.xlabel(r'$\\beta$'+'-coefficient in LASSO model')\n",
    "plt.ylabel('')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Both sex model\n",
    "tempDF = metWHtR_B_bcoefs\n",
    "bwhtr_color = 'b'\n",
    "bwhtr = 'MetWHtR'\n",
    "sex = 'BothSex'\n",
    "\n",
    "#Variables with non-zero beta-coefficient in all 10 models\n",
    "tempDF = tempDF.loc[tempDF['nZeros']==0].sort_values(by='Mean', ascending=False)\n",
    "tempDF = tempDF.drop(columns=['Mean', 'SD', 'nZeros'])\n",
    "print('Variables with non-zero beta-coefficient in all 10 models:', len(tempDF))\n",
    "\n",
    "#tidyr::gather in R\n",
    "tempDF1 = tempDF.reset_index().melt(var_name='LASSO', value_name='bcoef', id_vars=['Variable'])\n",
    "\n",
    "#Plot\n",
    "sns.set(style='ticks', font='Arial', context='talk')\n",
    "plt.figure(figsize=(5, 20))\n",
    "p = sns.boxplot(data=tempDF1, y='Variable', x='bcoef', color=bwhtr_color, dodge=False, saturation=1,\n",
    "                showfliers=True, flierprops={'marker':'o', 'markerfacecolor':'gray', 'alpha':0.4},\n",
    "                showcaps=True, notch=False)\n",
    "p.set(xlim=(-0.046, 0.046), xticks=np.arange(-0.04, 0.041, 0.02))#Fixed across omics\n",
    "p.grid(axis='x', linestyle='--', color='black')\n",
    "sns.despine()\n",
    "##Change default dull line color of sns.boxplot (saturation parameter is for patch)\n",
    "for line in p.get_lines():\n",
    "    line.set_color('k')\n",
    "for box in p.artists:\n",
    "    box.set_edgecolor('k')\n",
    "##Add background color\n",
    "for row_i in range(len(tempDF)):\n",
    "    if row_i%2 == 0:\n",
    "        plt.axhspan(ymin=row_i-0.5, ymax=row_i+0.5, facecolor=bwhtr_color, alpha=0.2, zorder=0)\n",
    "plt.ylabel('')\n",
    "plt.xlabel(r'$\\beta$'+'-coefficient in LASSO model\\n [log-scaled WHtR per s.d.]')\n",
    "##Save\n",
    "fileDir = './ExportFigures/'\n",
    "ipynbName = '220823_Multiomics-BMI-NatMed1stRevision_WHtR-LASSO-bcoef-ver2_'\n",
    "fileName = bwhtr+'-'+sex+'-bcoef_non-zero-in-all.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-5. Correlation b/w pairwise variables retained in all 10 models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Both sex model\n",
    "tempDF1 = metWHtR_B_bcoefs\n",
    "tempDF2 = metDF_B\n",
    "bwhtr_color = 'b'\n",
    "\n",
    "print('Both sex model:')\n",
    "#Variables with non-zero beta-coefficient in all 10 models\n",
    "tempL = tempDF1.loc[tempDF1['nZeros']==0].index.tolist()\n",
    "\n",
    "#Compute correlation matrix\n",
    "tempDF = tempDF2[tempL]\n",
    "tempDF = tempDF.corr(method='pearson')\n",
    "print('• Combinations:', int(len(tempDF)*(len(tempDF)-1)/2))\n",
    "\n",
    "#Extract lower triangle matrix\n",
    "tempDF = tempDF.where(np.tril(np.ones(tempDF.shape), k=-1).astype(np.bool), other=np.nan)\n",
    "tempDF.index.rename('Variable1', inplace=True)\n",
    "tempDF = tempDF.reset_index().melt(var_name='Variable2', value_name='Pearson_r', id_vars=['Variable1'])\n",
    "tempDF = tempDF.dropna()\n",
    "print('• nrows:', len(tempDF))\n",
    "print('• |Pearson\\'s r| > 0.8:', len(tempDF.loc[abs(tempDF['Pearson_r'])>0.8]))\n",
    "display(tempDF.loc[abs(tempDF['Pearson_r'])>0.8])\n",
    "\n",
    "#Distribution\n",
    "sns.set(style='ticks', font='Arial', context='talk')\n",
    "plt.figure(figsize=(4, 3))\n",
    "sns.distplot(tempDF['Pearson_r'], color=bwhtr_color).set(xlim=(-1, 1), xticks=np.arange(-1, 1.1, 0.5))\n",
    "plt.axvline(x=tempDF['Pearson_r'].max(), **{'linestyle':'--', 'color':bwhtr_color})\n",
    "plt.axvline(x=tempDF['Pearson_r'].min(), **{'linestyle':'--', 'color':bwhtr_color})\n",
    "sns.despine()\n",
    "plt.ylabel('Density')\n",
    "plt.xlabel('Pearson\\'s '+r'$r$')\n",
    "plt.show()\n",
    "\n",
    "print('Variables with non-zero beta-coefficient in all 10 models:', len(tempL))\n",
    "#Clustermap\n",
    "tempDF = tempDF2[tempL]\n",
    "tempDF = tempDF.corr(method='pearson')\n",
    "sns.set(style='ticks', font='Arial', context='notebook')\n",
    "cmap = sns.diverging_palette(220, 20, as_cmap=True)\n",
    "cm = sns.clustermap(tempDF, method='ward', metric='euclidean', cmap=cmap,\n",
    "                    row_cluster=True, col_cluster=True, row_linkage=None, col_linkage=None,\n",
    "                    row_colors=None, col_colors=None, xticklabels=True, yticklabels=True,\n",
    "                    dendrogram_ratio=(0.1, 0.1), cbar_pos=(0.9, 0.02, 0.02, 0.1),\n",
    "                    figsize=(15, 15), **{'center':0, 'vmin':-1, 'vmax':1})\n",
    "cm.cax.set_title('Pearson\\'s '+r'$r$')\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.5])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2-6. Explained variance in WHtR by the variables retained in ≥1 model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> OLS linear regression model for significance: log(WHtR) ~ b0 + b1\\*Variable + b2\\*Sex + b3\\*BaseAge + b4\\*PCs  \n",
    "> Univariate model: log(WHtR) ~ b0 + b1\\*Variable  \n",
    ">\n",
    "> Of note, biological WHtR is included for comparison.  \n",
    "> As well as LASSO models, unstandardized WHtR is used for the dependent variable.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempD1 = {'Female':metWHtR_F_bcoefs, 'Male':metWHtR_M_bcoefs, 'BothSex':metWHtR_B_bcoefs}\n",
    "tempD2 = {'Female':metDF_F, 'Male':metDF_M, 'BothSex':metDF_B}\n",
    "tempD3 = {'Female':whtrDF_F, 'Male':whtrDF_M, 'BothSex':whtrDF_B}\n",
    "yvar = 'WHtR'\n",
    "bwhtr = 'MetWHtR'\n",
    "topX = 30\n",
    "for sex in tempD1.keys():\n",
    "    print(sex+' model:')\n",
    "    \n",
    "    #Variables with non-zero beta-coefficient at least 1 of 10 models\n",
    "    tempDF = tempD1[sex]\n",
    "    tempL = tempDF.loc[tempDF['nZeros']==10].index.tolist()#To be excluded\n",
    "    \n",
    "    #Prepare DF for regressions\n",
    "    tempDF = tempD2[sex].drop(columns=tempL)\n",
    "    tempDF = tempDF.rename(columns={'log_Base'+bwhtr:bwhtr+' model'})\n",
    "    #Replace the standardized dependent variable with the unstandardized\n",
    "    tempDF1 = tempD3[sex]\n",
    "    tempDF['log_Base'+yvar] = tempDF1['log_Base'+yvar]\n",
    "    \n",
    "    #Independent variables for regressions (including bWHtR)\n",
    "    tempL = tempDF.drop(columns=whtrDF.columns.tolist()).columns.tolist()\n",
    "    \n",
    "    #Perform OLS linear regression\n",
    "    tempL1 = []#For sample size\n",
    "    tempL2 = []#For R2 in the univariate model\n",
    "    tempL3 = []#For R2\n",
    "    tempL4 = []#For beta-coefficient\n",
    "    tempL5 = []#For SE of beta-coefficient\n",
    "    tempL6 = []#For t-statistic\n",
    "    tempL7 = []#For residual degrees of freedom\n",
    "    tempL8 = []#For P-value\n",
    "    t_start = time.time()\n",
    "    for var in tempL:\n",
    "        #Rename independent variable\n",
    "        tempDF1 = tempDF.rename(columns={var: 'Variable'})\n",
    "        #Add a constant for the intercept -> Similar to R, smf automatically add a constant\n",
    "        \n",
    "        #Save sample size\n",
    "        tempL1.append(len(tempDF1))\n",
    "        \n",
    "        #Fit univariate model\n",
    "        formula = 'log_Base'+yvar+' ~ Variable'\n",
    "        fit_res = smf.ols(formula, data=tempDF1).fit()\n",
    "        #Save R2 [%]\n",
    "        tempL2.append(fit_res.rsquared*100)\n",
    "\n",
    "        #Fit full model\n",
    "        if sex=='BothSex':\n",
    "            formula = 'log_Base'+yvar+' ~ Variable + C(Sex) + BaseAge + PC1 + PC2 + PC3 + PC4 + PC5'\n",
    "        else:#Sex is invariant\n",
    "            formula = 'log_Base'+yvar+' ~ Variable + BaseAge + PC1 + PC2 + PC3 + PC4 + PC5'\n",
    "        fit_res = smf.ols(formula, data=tempDF1).fit()\n",
    "        #Save R2 [%]\n",
    "        tempL3.append(fit_res.rsquared*100)\n",
    "        #Save beta-coefficient of the variable\n",
    "        tempL4.append(fit_res.params['Variable'])\n",
    "        tempL5.append(fit_res.bse['Variable'])\n",
    "        #Save t-statistic of the variable\n",
    "        tempL6.append(fit_res.tvalues['Variable'])\n",
    "        #Save residual degrees of freedom\n",
    "        tempL7.append(int(fit_res.df_resid))\n",
    "        #Save P-value of the variable\n",
    "        tempL8.append(fit_res.pvalues['Variable'])\n",
    "    t_elapsed = time.time() - t_start\n",
    "    print('Elapsed time for', len(tempL), 'OLS linear regressions (including '+bwhtr+'):',\n",
    "          round(t_elapsed//60), 'min', round(t_elapsed%60, 1), 'sec')\n",
    "    \n",
    "    #Clean the results\n",
    "    tempDF = pd.DataFrame({'N':tempL1, 'UnivarR2':tempL2, 'R2':tempL3,\n",
    "                           'OLSbcoef':tempL4, 'OLSbcoefSE':tempL5,\n",
    "                           'tStat':tempL6, 'DoF':tempL7,  'Pval':tempL8},\n",
    "                          index=pd.Index(tempL, name='Variable'))\n",
    "    ##P-value adjustment by using Benjamini–Hochberg method\n",
    "    tempDF['AdjPval'] = multi.multipletests(tempDF['Pval'], alpha=0.05, method='fdr_bh',\n",
    "                                            is_sorted=False, returnsorted=False)[1]\n",
    "    ##Add the LASSO results\n",
    "    tempDF1 = tempD1[sex]\n",
    "    tempDF['LASSObcoef'] = tempDF1.loc[tempDF1.index.isin(tempL)]['Mean']#NaN in bWHtR reference\n",
    "    tempDF['LASSOnZeros'] = tempDF1.loc[tempDF1.index.isin(tempL)]['nZeros']#NaN in bWHtR reference\n",
    "    tempDF = tempDF.sort_values(by='AdjPval', ascending=True)\n",
    "    \n",
    "    #Save the cleaned DF\n",
    "    fileDir = './ExportData/'\n",
    "    ipynbName = '220823_Multiomics-BMI-NatMed1stRevision_WHtR-LASSO-bcoef-ver2_'\n",
    "    fileName = bwhtr+'-'+sex+'-OLS.tsv'\n",
    "    tempDF.to_csv(fileDir+ipynbName+fileName, sep='\\t', index=True)\n",
    "    \n",
    "    #Extact significant variables\n",
    "    tempDF = tempDF[tempDF['AdjPval']<0.05]\n",
    "    print('Variables significantly associated with WHtR (FDR < 0.05):', len(tempDF)-1)#bWHtR reference\n",
    "    ##Top X (+ bWHtR)\n",
    "    tempDF = tempDF.iloc[0:(topX+1)]\n",
    "    display(tempDF)\n",
    "    ##Category and color\n",
    "    tempL = []\n",
    "    for row_n in tempDF.index.tolist():\n",
    "        if row_n==bwhtr+' model':\n",
    "            tempL.append('Positive association')\n",
    "        else:\n",
    "            bcoef = tempDF.loc[row_n, 'OLSbcoef']\n",
    "            if bcoef>0:\n",
    "                tempL.append('Positive association')\n",
    "            elif bcoef<0:\n",
    "                tempL.append('Negative association')\n",
    "            else:\n",
    "                tempL.append('Error?')\n",
    "    tempDF['Association'] = tempL\n",
    "    tempD = {'Positive association':'tab:red', 'Negative association':'tab:blue'}\n",
    "    ##Plot\n",
    "    tempDF = tempDF.reset_index()\n",
    "    sns.set(style='ticks', font='Arial', context='talk')\n",
    "    plt.figure(figsize=(10, 4))\n",
    "    p = sns.barplot(data=tempDF, y='UnivarR2', x='Variable', hue='Association', hue_order=tempD.keys(),\n",
    "                    dodge=False, palette=tempD, edgecolor='k')\n",
    "    p.set(ylim=(0, 80), yticks=np.arange(0, 80, 15))#Fixed across omics\n",
    "    p.grid(axis='y', linestyle='--', color='black')\n",
    "    sns.despine()\n",
    "    plt.xticks(rotation=90, horizontalalignment='center')\n",
    "    plt.xlabel('')\n",
    "    plt.ylabel('Explained variance in WHtR [%]')\n",
    "    plt.legend(loc='upper right')\n",
    "    if sex=='BothSex':#Save\n",
    "        fileDir = './ExportFigures/'\n",
    "        ipynbName = '220823_Multiomics-BMI-NatMed1stRevision_WHtR-LASSO-bcoef-ver2_'\n",
    "        fileName = bwhtr+'-'+sex+'-OLS-top'+str(topX)+'.tif'\n",
    "        plt.gcf().savefig(fileDir+ipynbName+fileName, dpi=300, bbox_inches='tight', pad_inches=0.04,\n",
    "                          pil_kwargs={'compression':'tiff_lzw'})\n",
    "    plt.show()\n",
    "    print('')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Proteomics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3-1. Prepare and standardize variables"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Variables, including bWHtR (as a reference) and covariates, are standardized for OLS linear regression.  \n",
    "> –> In this version, unstandardized WHtR is used as the dependent variable. However, WHtR within the following DFs is standardized (just for the code simplicity).  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempDF = protDF\n",
    "tempD1 = {'Female':whtrDF_F, 'Male':whtrDF_M, 'BothSex':whtrDF_B}\n",
    "bwhtr = 'ProtWHtR'\n",
    "tempD2 = {}\n",
    "for sex in tempD1.keys():\n",
    "    #Add covariates (and WHtR) while selecting the cohort\n",
    "    tempDF1 = tempD1[sex]\n",
    "    tempDF1 = pd.merge(tempDF, tempDF1, left_index=True, right_index=True, how='inner')\n",
    "    \n",
    "    #Add biological WHtR\n",
    "    fileDir = './ExportData/'\n",
    "    ipynbName = '220822_Multiomics-BMI-NatMed1stRevision_WHtR-baseline-LASSO-ver2_'\n",
    "    fileName = bwhtr+'-'+sex+'.tsv'\n",
    "    tempDF2 = pd.read_csv(fileDir+ipynbName+fileName, sep='\\t', dtype={'public_client_id':str})\n",
    "    tempDF2 = tempDF2.set_index('public_client_id')\n",
    "    tempDF1['log_Base'+bwhtr] = tempDF2['log_Base'+bwhtr]\n",
    "    \n",
    "    #Z-score transformation\n",
    "    tempDF2 = tempDF1.select_dtypes(include=[np.number])\n",
    "    scaler = StandardScaler(copy=True, with_mean=True, with_std=True)\n",
    "    tempA = scaler.fit_transform(tempDF2)#Column direction\n",
    "    tempDF2 = pd.DataFrame(data=tempA, index=tempDF2.index, columns=tempDF2.columns)\n",
    "    ##Recover categorical covariates\n",
    "    tempDF3 = tempDF1.select_dtypes(exclude=[np.number])\n",
    "    tempDF1 = pd.merge(tempDF2, tempDF3, left_index=True, right_index=True, how='left')\n",
    "    \n",
    "    print(sex)\n",
    "    display(tempDF1.describe(include='all'))\n",
    "    tempD2[sex] = tempDF1\n",
    "\n",
    "protDF_F = tempD2['Female']\n",
    "protDF_M = tempD2['Male']\n",
    "protDF_B = tempD2['BothSex']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3-2. Correlation b/w all pairwise variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Female and male cohorts\n",
    "tempD1 = {'Female':protDF_F, 'Male':protDF_M}\n",
    "bwhtr = 'ProtWHtR'\n",
    "tempD2 = {}\n",
    "for sex in tempD1.keys():\n",
    "    print(sex+' cohort:')\n",
    "    tempDF = tempD1[sex]\n",
    "    #Compute correlation matrix\n",
    "    tempL = [item for sublist in [whtrDF.columns.tolist(), ['log_Base'+bwhtr]] for item in sublist]\n",
    "    tempDF1 = tempDF.drop(columns=tempL)\n",
    "    tempDF1 = tempDF1.corr(method='pearson')\n",
    "    print('• Combinations:', int(len(tempDF1)*(len(tempDF1)-1)/2))\n",
    "    \n",
    "    #Extract lower triangle matrix\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",
    "    print('• nrows:', len(tempDF1))\n",
    "    print('• |Pearson\\'s r| > 0.8:', len(tempDF1.loc[abs(tempDF1['Pearson_r'])>0.8]))\n",
    "    display(tempDF1.loc[abs(tempDF1['Pearson_r'])>0.8])\n",
    "    \n",
    "    tempD2[sex] = tempDF1\n",
    "\n",
    "#Distribution\n",
    "tempD3 = {'Female':'tab:red', 'Male':'tab:blue'}\n",
    "sns.set(style='ticks', font='Arial', context='notebook')\n",
    "plt.figure(figsize=(4, 3))\n",
    "for sex in tempD2.keys():\n",
    "    tempS = tempD2[sex]['Pearson_r']\n",
    "    sns.distplot(tempS, color=tempD3[sex], label=sex).set(xlim=(-1, 1), xticks=np.arange(-1, 1.1, 0.5))\n",
    "sns.despine()\n",
    "plt.xlabel('Pearson\\'s '+r'$r$')\n",
    "plt.ylabel('Density')\n",
    "plt.legend(loc='upper right')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Both sex cohort\n",
    "tempDF = protDF_B\n",
    "bwhtr = 'ProtWHtR'\n",
    "bwhtr_color = 'r'\n",
    "\n",
    "print('Both sex cohort:')\n",
    "#Compute correlation matrix\n",
    "tempL = [item for sublist in [whtrDF.columns.tolist(), ['log_Base'+bwhtr]] for item in sublist]\n",
    "tempDF = tempDF.drop(columns=tempL)\n",
    "tempDF = tempDF.corr(method='pearson')\n",
    "print('• Combinations:', int(len(tempDF)*(len(tempDF)-1)/2))\n",
    "\n",
    "#Extract lower triangle matrix\n",
    "tempDF = tempDF.where(np.tril(np.ones(tempDF.shape), k=-1).astype(np.bool), other=np.nan)\n",
    "tempDF.index.rename('Variable1', inplace=True)\n",
    "tempDF = tempDF.reset_index().melt(var_name='Variable2', value_name='Pearson_r', id_vars=['Variable1'])\n",
    "tempDF = tempDF.dropna()\n",
    "print('• nrows:', len(tempDF))\n",
    "print('• |Pearson\\'s r| > 0.8:', len(tempDF.loc[abs(tempDF['Pearson_r'])>0.8]))\n",
    "display(tempDF.loc[abs(tempDF['Pearson_r'])>0.8])\n",
    "\n",
    "#Distribution\n",
    "sns.set(style='ticks', font='Arial', context='talk')\n",
    "plt.figure(figsize=(4, 3))\n",
    "sns.distplot(tempDF['Pearson_r'], color=bwhtr_color).set(xlim=(-1, 1), xticks=np.arange(-1, 1.1, 0.5))\n",
    "plt.axvline(x=tempDF['Pearson_r'].max(), **{'linestyle':'--', 'color':bwhtr_color})\n",
    "plt.axvline(x=tempDF['Pearson_r'].min(), **{'linestyle':'--', 'color':bwhtr_color})\n",
    "sns.despine()\n",
    "plt.ylabel('Density')\n",
    "plt.xlabel('Pearson\\'s '+r'$r$')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3-3. Prepare the beta-coefficients of LASSO models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bwhtr = 'ProtWHtR'\n",
    "tempD = {}\n",
    "for sex in ['Female', 'Male', 'BothSex']:\n",
    "    #Import the LASSO beta-coefficients\n",
    "    fileDir = './ExportData/'\n",
    "    ipynbName = '220822_Multiomics-BMI-NatMed1stRevision_WHtR-baseline-LASSO-ver2_'\n",
    "    fileName = bwhtr+'-'+sex+'-LASSObcoefs.tsv'\n",
    "    tempDF = pd.read_csv(fileDir+ipynbName+fileName, sep='\\t').set_index('Variable')\n",
    "    tempDF = tempDF.drop(index=['Intercept'])\n",
    "    tempD[sex] = tempDF\n",
    "    \n",
    "    #Check\n",
    "    print(sex+':')\n",
    "    print(' - Variables:', len(tempDF))\n",
    "    tempDF1 = tempDF.loc[tempDF['nZeros']!=10]\n",
    "    print(' - Variables with non-zero beta-coefficient:', len(tempDF1),\n",
    "          '(', len(tempDF1)/len(tempDF)*100, '%)')\n",
    "    tempDF1 = tempDF.loc[tempDF['nZeros']==0]\n",
    "    print(' - Variables with non-zero beta-coefficient in all 10 models:', len(tempDF1),\n",
    "          '(', len(tempDF1)/len(tempDF)*100, '%)')\n",
    "    tempDF1 = tempDF1.sort_values(by='Mean', ascending=False)\n",
    "    display(tempDF1)\n",
    "    print('')\n",
    "\n",
    "protWHtR_F_bcoefs = tempD['Female']\n",
    "protWHtR_M_bcoefs = tempD['Male']\n",
    "protWHtR_B_bcoefs = tempD['BothSex']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3-4. Variables retained in all 10 models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempDF1 = protWHtR_F_bcoefs\n",
    "tempDF2 = protWHtR_M_bcoefs\n",
    "tempDF3 = protWHtR_B_bcoefs\n",
    "category = 'Proteomics'\n",
    "\n",
    "#Variables with non-zero beta-coefficient in all 10 models\n",
    "tempS1 = set(tempDF1.loc[tempDF1['nZeros']==0].index)\n",
    "tempS2 = set(tempDF2.loc[tempDF2['nZeros']==0].index)\n",
    "tempS3 = set(tempDF3.loc[tempDF3['nZeros']==0].index)\n",
    "\n",
    "#Common the robust variables\n",
    "var_111 = tempS1 & tempS2 & tempS3\n",
    "print('Common variables with non-zero beta-coefficient in all 10 models:\\n', var_111)\n",
    "var_110 = (tempS1 & tempS2) - var_111\n",
    "var_101 = (tempS1 & tempS3) - var_111\n",
    "var_011 = (tempS2 & tempS3) - var_111\n",
    "var_100 = tempS1 - (var_110 | var_101 | var_111)\n",
    "var_010 = tempS2 - (var_110 | var_011 | var_111)\n",
    "var_001 = tempS3 - (var_101 | var_011 | var_111)\n",
    "\n",
    "#The number of each subset\n",
    "var_subset = [len(var_100), len(var_010), len(var_110), len(var_001),\n",
    "              len(var_101), len(var_011), len(var_111)]\n",
    "print('Each subset:', var_subset)\n",
    "\n",
    "#Venn diagram\n",
    "sns.set(font='Arial', context='talk')\n",
    "venn3(subsets=var_subset, set_labels=('Female', 'Male', 'Both sex'),\n",
    "      set_colors=('tab:red', 'tab:blue', 'tab:green'), alpha=0.4)\n",
    "venn3_circles(subsets=var_subset)\n",
    "plt.title('Robust variables\\n—'+category+'—', fontdict={'fontsize':24})\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempDF1 = protWHtR_F_bcoefs\n",
    "tempDF2 = protWHtR_M_bcoefs\n",
    "tempDF3 = protWHtR_B_bcoefs\n",
    "bwhtr_color = 'r'\n",
    "\n",
    "#Variables with non-zero beta-coefficient in all 10 models\n",
    "tempS1 = set(tempDF1.loc[tempDF1['nZeros']==0].index)\n",
    "tempS2 = set(tempDF2.loc[tempDF2['nZeros']==0].index)\n",
    "tempS3 = set(tempDF3.loc[tempDF3['nZeros']==0].index)\n",
    "\n",
    "#Union of the robust variables\n",
    "tempL = list(tempS1 | tempS2 | tempS3)\n",
    "print('Union of the variables with non-zero beta-coefficient in all 10 models:', len(tempL))\n",
    "\n",
    "#Merge (use original bcoefs to rescue ones with zeros in any 10 models)\n",
    "tempDF1 = tempDF1.loc[tempDF1.index.isin(tempL)].drop(columns=['Mean', 'SD', 'nZeros'])\n",
    "tempDF2 = tempDF2.loc[tempDF2.index.isin(tempL)].drop(columns=['Mean', 'SD', 'nZeros'])\n",
    "tempDF3 = tempDF3.loc[tempDF3.index.isin(tempL)].drop(columns=['Mean', 'SD', 'nZeros'])\n",
    "tempDF = pd.concat([tempDF1, tempDF2, tempDF3], axis=0)\n",
    "tempDF['Model'] = np.repeat(['Female', 'Male', 'Both sex'], [len(tempDF1), len(tempDF2), len(tempDF3)])\n",
    "\n",
    "#tidyr::gather in R\n",
    "tempDF1 = tempDF.reset_index().melt(var_name='LASSO', value_name='bcoef', id_vars=['Variable', 'Model'])\n",
    "\n",
    "#Order for plot\n",
    "tempDF2 = tempDF1.loc[tempDF1['Model']=='Both sex'].groupby('Variable', as_index=True).agg({'bcoef':np.mean})\n",
    "tempDF2 = tempDF2.sort_values(by='bcoef', ascending=False)\n",
    "\n",
    "#Plot\n",
    "sns.set(style='ticks', font='Arial', context='talk')\n",
    "sns.catplot(data=tempDF1, y='Variable', order=tempDF2.index.tolist(), x='bcoef', kind='box',\n",
    "            hue='Model', palette='Set1', dodge=True, height=14, aspect=0.6,\n",
    "            showfliers=True, flierprops={'marker':'o', 'markerfacecolor':'gray', 'alpha':0.4})\n",
    "plt.grid(axis='x', linestyle='--', color='black')\n",
    "for row_i in range(len(tempDF2)):\n",
    "    if row_i%2 == 0:\n",
    "        plt.axhspan(ymin=row_i-0.5, ymax=row_i+0.5, facecolor=bwhtr_color, alpha=0.2, zorder=0)\n",
    "plt.xlabel(r'$\\beta$'+'-coefficient in LASSO model')\n",
    "plt.ylabel('')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Both sex model\n",
    "tempDF = protWHtR_B_bcoefs\n",
    "bwhtr_color = 'r'\n",
    "bwhtr = 'ProtWHtR'\n",
    "sex = 'BothSex'\n",
    "\n",
    "#Variables with non-zero beta-coefficient in all 10 models\n",
    "tempDF = tempDF.loc[tempDF['nZeros']==0].sort_values(by='Mean', ascending=False)\n",
    "tempDF = tempDF.drop(columns=['Mean', 'SD', 'nZeros'])\n",
    "print('Variables with non-zero beta-coefficient in all 10 models:', len(tempDF))\n",
    "\n",
    "#tidyr::gather in R\n",
    "tempDF1 = tempDF.reset_index().melt(var_name='LASSO', value_name='bcoef', id_vars=['Variable'])\n",
    "\n",
    "#Plot\n",
    "sns.set(style='ticks', font='Arial', context='talk')\n",
    "plt.figure(figsize=(5, 9))\n",
    "p = sns.boxplot(data=tempDF1, y='Variable', x='bcoef', color=bwhtr_color, dodge=False, saturation=1,\n",
    "                showfliers=True, flierprops={'marker':'o', 'markerfacecolor':'gray', 'alpha':0.4},\n",
    "                showcaps=True, notch=False)\n",
    "p.set(xlim=(-0.046, 0.046), xticks=np.arange(-0.04, 0.041, 0.02))#Fixed across omics\n",
    "p.grid(axis='x', linestyle='--', color='black')\n",
    "sns.despine()\n",
    "##Change default dull line color of sns.boxplot (saturation parameter is for patch)\n",
    "for line in p.get_lines():\n",
    "    line.set_color('k')\n",
    "for box in p.artists:\n",
    "    box.set_edgecolor('k')\n",
    "##Add background color\n",
    "for row_i in range(len(tempDF)):\n",
    "    if row_i%2 == 0:\n",
    "        plt.axhspan(ymin=row_i-0.5, ymax=row_i+0.5, facecolor=bwhtr_color, alpha=0.2, zorder=0)\n",
    "plt.ylabel('')\n",
    "plt.xlabel(r'$\\beta$'+'-coefficient in LASSO model\\n [log-scaled WHtR per s.d.]')\n",
    "##Save\n",
    "fileDir = './ExportFigures/'\n",
    "ipynbName = '220823_Multiomics-BMI-NatMed1stRevision_WHtR-LASSO-bcoef-ver2_'\n",
    "fileName = bwhtr+'-'+sex+'-bcoef_non-zero-in-all.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": [
    "### 3-5. Correlation b/w pairwise variables retained in all 10 models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Both sex model\n",
    "tempDF1 = protWHtR_B_bcoefs\n",
    "tempDF2 = protDF_B\n",
    "bwhtr_color = 'r'\n",
    "\n",
    "print('Both sex model:')\n",
    "#Variables with non-zero beta-coefficient in all 10 models\n",
    "tempL = tempDF1.loc[tempDF1['nZeros']==0].index.tolist()\n",
    "\n",
    "#Compute correlation matrix\n",
    "tempDF = tempDF2[tempL]\n",
    "tempDF = tempDF.corr(method='pearson')\n",
    "print('• Combinations:', int(len(tempDF)*(len(tempDF)-1)/2))\n",
    "\n",
    "#Extract lower triangle matrix\n",
    "tempDF = tempDF.where(np.tril(np.ones(tempDF.shape), k=-1).astype(np.bool), other=np.nan)\n",
    "tempDF.index.rename('Variable1', inplace=True)\n",
    "tempDF = tempDF.reset_index().melt(var_name='Variable2', value_name='Pearson_r', id_vars=['Variable1'])\n",
    "tempDF = tempDF.dropna()\n",
    "print('• nrows:', len(tempDF))\n",
    "print('• |Pearson\\'s r| > 0.8:', len(tempDF.loc[abs(tempDF['Pearson_r'])>0.8]))\n",
    "display(tempDF.loc[abs(tempDF['Pearson_r'])>0.8])\n",
    "\n",
    "#Distribution\n",
    "sns.set(style='ticks', font='Arial', context='talk')\n",
    "plt.figure(figsize=(4, 3))\n",
    "sns.distplot(tempDF['Pearson_r'], color=bwhtr_color).set(xlim=(-1, 1), xticks=np.arange(-1, 1.1, 0.5))\n",
    "plt.axvline(x=tempDF['Pearson_r'].max(), **{'linestyle':'--', 'color':bwhtr_color})\n",
    "plt.axvline(x=tempDF['Pearson_r'].min(), **{'linestyle':'--', 'color':bwhtr_color})\n",
    "sns.despine()\n",
    "plt.ylabel('Density')\n",
    "plt.xlabel('Pearson\\'s '+r'$r$')\n",
    "plt.show()\n",
    "\n",
    "print('Variables with non-zero beta-coefficient in all 10 models:', len(tempL))\n",
    "#Clustermap\n",
    "tempDF = tempDF2[tempL]\n",
    "tempDF = tempDF.corr(method='pearson')\n",
    "sns.set(style='ticks', font='Arial', context='notebook')\n",
    "cmap = sns.diverging_palette(220, 20, as_cmap=True)\n",
    "cm = sns.clustermap(tempDF, method='ward', metric='euclidean', cmap=cmap,\n",
    "                    row_cluster=True, col_cluster=True, row_linkage=None, col_linkage=None,\n",
    "                    row_colors=None, col_colors=None, xticklabels=True, yticklabels=True,\n",
    "                    dendrogram_ratio=(0.15, 0.15), cbar_pos=(0.91, 0.025, 0.02, 0.06),\n",
    "                    figsize=(7, 7), **{'center':0, 'vmin':-1, 'vmax':1})\n",
    "cm.cax.set_title('Pearson\\'s '+r'$r$')\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.5])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3-6. Explained variance in WHtR by the variables retained in ≥1 model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> OLS linear regression model for significance: log(WHtR) ~ b0 + b1\\*Variable + b2\\*Sex + b3\\*BaseAge + b4\\*PCs  \n",
    "> Univariate model: log(WHtR) ~ b0 + b1\\*Variable  \n",
    ">\n",
    "> Of note, biological WHtR is included for comparison.  \n",
    "> As well as LASSO models, unstandardized WHtR is used for the dependent variable.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempD1 = {'Female':protWHtR_F_bcoefs, 'Male':protWHtR_M_bcoefs, 'BothSex':protWHtR_B_bcoefs}\n",
    "tempD2 = {'Female':protDF_F, 'Male':protDF_M, 'BothSex':protDF_B}\n",
    "tempD3 = {'Female':whtrDF_F, 'Male':whtrDF_M, 'BothSex':whtrDF_B}\n",
    "yvar = 'WHtR'\n",
    "bwhtr = 'ProtWHtR'\n",
    "topX = 30\n",
    "for sex in tempD1.keys():\n",
    "    print(sex+' model:')\n",
    "    \n",
    "    #Variables with non-zero beta-coefficient at least 1 of 10 models\n",
    "    tempDF = tempD1[sex]\n",
    "    tempL = tempDF.loc[tempDF['nZeros']==10].index.tolist()#To be excluded\n",
    "    \n",
    "    #Prepare DF for regressions\n",
    "    tempDF = tempD2[sex].drop(columns=tempL)\n",
    "    tempDF = tempDF.rename(columns={'log_Base'+bwhtr:bwhtr+' model'})\n",
    "    #Replace the standardized dependent variable with the unstandardized\n",
    "    tempDF1 = tempD3[sex]\n",
    "    tempDF['log_Base'+yvar] = tempDF1['log_Base'+yvar]\n",
    "    \n",
    "    #Independent variables for regressions (including bWHtR)\n",
    "    tempL = tempDF.drop(columns=whtrDF.columns.tolist()).columns.tolist()\n",
    "    \n",
    "    #Perform OLS linear regression\n",
    "    tempL1 = []#For sample size\n",
    "    tempL2 = []#For R2 in the univariate model\n",
    "    tempL3 = []#For R2\n",
    "    tempL4 = []#For beta-coefficient\n",
    "    tempL5 = []#For SE of beta-coefficient\n",
    "    tempL6 = []#For t-statistic\n",
    "    tempL7 = []#For residual degrees of freedom\n",
    "    tempL8 = []#For P-value\n",
    "    t_start = time.time()\n",
    "    for var in tempL:\n",
    "        #Rename independent variable\n",
    "        tempDF1 = tempDF.rename(columns={var: 'Variable'})\n",
    "        #Add a constant for the intercept -> Similar to R, smf automatically add a constant\n",
    "        \n",
    "        #Save sample size\n",
    "        tempL1.append(len(tempDF1))\n",
    "        \n",
    "        #Fit univariate model\n",
    "        formula = 'log_Base'+yvar+' ~ Variable'\n",
    "        fit_res = smf.ols(formula, data=tempDF1).fit()\n",
    "        #Save R2 [%]\n",
    "        tempL2.append(fit_res.rsquared*100)\n",
    "\n",
    "        #Fit full model\n",
    "        if sex=='BothSex':\n",
    "            formula = 'log_Base'+yvar+' ~ Variable + C(Sex) + BaseAge + PC1 + PC2 + PC3 + PC4 + PC5'\n",
    "        else:#Sex is invariant\n",
    "            formula = 'log_Base'+yvar+' ~ Variable + BaseAge + PC1 + PC2 + PC3 + PC4 + PC5'\n",
    "        fit_res = smf.ols(formula, data=tempDF1).fit()\n",
    "        #Save R2 [%]\n",
    "        tempL3.append(fit_res.rsquared*100)\n",
    "        #Save beta-coefficient of the variable\n",
    "        tempL4.append(fit_res.params['Variable'])\n",
    "        tempL5.append(fit_res.bse['Variable'])\n",
    "        #Save t-statistic of the variable\n",
    "        tempL6.append(fit_res.tvalues['Variable'])\n",
    "        #Save residual degrees of freedom\n",
    "        tempL7.append(int(fit_res.df_resid))\n",
    "        #Save P-value of the variable\n",
    "        tempL8.append(fit_res.pvalues['Variable'])\n",
    "    t_elapsed = time.time() - t_start\n",
    "    print('Elapsed time for', len(tempL), 'OLS linear regressions (including '+bwhtr+'):',\n",
    "          round(t_elapsed//60), 'min', round(t_elapsed%60, 1), 'sec')\n",
    "    \n",
    "    #Clean the results\n",
    "    tempDF = pd.DataFrame({'N':tempL1, 'UnivarR2':tempL2, 'R2':tempL3,\n",
    "                           'OLSbcoef':tempL4, 'OLSbcoefSE':tempL5,\n",
    "                           'tStat':tempL6, 'DoF':tempL7,  'Pval':tempL8},\n",
    "                          index=pd.Index(tempL, name='Variable'))\n",
    "    ##P-value adjustment by using Benjamini–Hochberg method\n",
    "    tempDF['AdjPval'] = multi.multipletests(tempDF['Pval'], alpha=0.05, method='fdr_bh',\n",
    "                                            is_sorted=False, returnsorted=False)[1]\n",
    "    ##Add the LASSO results\n",
    "    tempDF1 = tempD1[sex]\n",
    "    tempDF['LASSObcoef'] = tempDF1.loc[tempDF1.index.isin(tempL)]['Mean']#NaN in bWHtR reference\n",
    "    tempDF['LASSOnZeros'] = tempDF1.loc[tempDF1.index.isin(tempL)]['nZeros']#NaN in bWHtR reference\n",
    "    tempDF = tempDF.sort_values(by='AdjPval', ascending=True)\n",
    "    \n",
    "    #Save the cleaned DF\n",
    "    fileDir = './ExportData/'\n",
    "    ipynbName = '220823_Multiomics-BMI-NatMed1stRevision_WHtR-LASSO-bcoef-ver2_'\n",
    "    fileName = bwhtr+'-'+sex+'-OLS.tsv'\n",
    "    tempDF.to_csv(fileDir+ipynbName+fileName, sep='\\t', index=True)\n",
    "    \n",
    "    #Extact significant variables\n",
    "    tempDF = tempDF[tempDF['AdjPval']<0.05]\n",
    "    print('Variables significantly associated with WHtR (FDR < 0.05):', len(tempDF)-1)#bWHtR reference\n",
    "    ##Top X (+ bWHtR)\n",
    "    tempDF = tempDF.iloc[0:(topX+1)]\n",
    "    display(tempDF)\n",
    "    ##Category and color\n",
    "    tempL = []\n",
    "    for row_n in tempDF.index.tolist():\n",
    "        if row_n==bwhtr+' model':\n",
    "            tempL.append('Positive association')\n",
    "        else:\n",
    "            bcoef = tempDF.loc[row_n, 'OLSbcoef']\n",
    "            if bcoef>0:\n",
    "                tempL.append('Positive association')\n",
    "            elif bcoef<0:\n",
    "                tempL.append('Negative association')\n",
    "            else:\n",
    "                tempL.append('Error?')\n",
    "    tempDF['Association'] = tempL\n",
    "    tempD = {'Positive association':'tab:red', 'Negative association':'tab:blue'}\n",
    "    ##Plot\n",
    "    tempDF = tempDF.reset_index()\n",
    "    sns.set(style='ticks', font='Arial', context='talk')\n",
    "    plt.figure(figsize=(10, 4))\n",
    "    p = sns.barplot(data=tempDF, y='UnivarR2', x='Variable', hue='Association', hue_order=tempD.keys(),\n",
    "                    dodge=False, palette=tempD, edgecolor='k')\n",
    "    p.set(ylim=(0, 80), yticks=np.arange(0, 80, 15))#Fixed across omics\n",
    "    p.grid(axis='y', linestyle='--', color='black')\n",
    "    sns.despine()\n",
    "    plt.xticks(rotation=90, horizontalalignment='center')\n",
    "    plt.xlabel('')\n",
    "    plt.ylabel('Explained variance in WHtR [%]')\n",
    "    plt.legend(loc='upper right')\n",
    "    if sex=='BothSex':#Save\n",
    "        fileDir = './ExportFigures/'\n",
    "        ipynbName = '220823_Multiomics-BMI-NatMed1stRevision_WHtR-LASSO-bcoef-ver2_'\n",
    "        fileName = bwhtr+'-'+sex+'-OLS-top'+str(topX)+'.tif'\n",
    "        plt.gcf().savefig(fileDir+ipynbName+fileName, dpi=300, bbox_inches='tight', pad_inches=0.04,\n",
    "                          pil_kwargs={'compression':'tiff_lzw'})\n",
    "    plt.show()\n",
    "    print('')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Clinical labs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4-1. Prepare and standardize variables"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Variables, including bWHtR (as a reference) and covariates, are standardized for OLS linear regression.  \n",
    "> –> In this version, unstandardized WHtR is used as the dependent variable. However, WHtR within the following DFs is standardized (just for the code simplicity).  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempDF = chemDF\n",
    "tempD1 = {'Female':whtrDF_F, 'Male':whtrDF_M, 'BothSex':whtrDF_B}\n",
    "bwhtr = 'ChemWHtR'\n",
    "tempD2 = {}\n",
    "for sex in tempD1.keys():\n",
    "    #Add covariates (and WHtR) while selecting the cohort\n",
    "    tempDF1 = tempD1[sex]\n",
    "    tempDF1 = pd.merge(tempDF, tempDF1, left_index=True, right_index=True, how='inner')\n",
    "    \n",
    "    #Add biological WHtR\n",
    "    fileDir = './ExportData/'\n",
    "    ipynbName = '220822_Multiomics-BMI-NatMed1stRevision_WHtR-baseline-LASSO-ver2_'\n",
    "    fileName = bwhtr+'-'+sex+'.tsv'\n",
    "    tempDF2 = pd.read_csv(fileDir+ipynbName+fileName, sep='\\t', dtype={'public_client_id':str})\n",
    "    tempDF2 = tempDF2.set_index('public_client_id')\n",
    "    tempDF1['log_Base'+bwhtr] = tempDF2['log_Base'+bwhtr]\n",
    "    \n",
    "    #Z-score transformation\n",
    "    tempDF2 = tempDF1.select_dtypes(include=[np.number])\n",
    "    scaler = StandardScaler(copy=True, with_mean=True, with_std=True)\n",
    "    tempA = scaler.fit_transform(tempDF2)#Column direction\n",
    "    tempDF2 = pd.DataFrame(data=tempA, index=tempDF2.index, columns=tempDF2.columns)\n",
    "    ##Recover categorical covariates\n",
    "    tempDF3 = tempDF1.select_dtypes(exclude=[np.number])\n",
    "    tempDF1 = pd.merge(tempDF2, tempDF3, left_index=True, right_index=True, how='left')\n",
    "    \n",
    "    print(sex)\n",
    "    display(tempDF1.describe(include='all'))\n",
    "    tempD2[sex] = tempDF1\n",
    "\n",
    "chemDF_F = tempD2['Female']\n",
    "chemDF_M = tempD2['Male']\n",
    "chemDF_B = tempD2['BothSex']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4-2. Correlation b/w all pairwise variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Female and male cohorts\n",
    "tempD1 = {'Female':chemDF_F, 'Male':chemDF_M}\n",
    "bwhtr = 'ChemWHtR'\n",
    "tempD2 = {}\n",
    "for sex in tempD1.keys():\n",
    "    print(sex+' cohort:')\n",
    "    tempDF = tempD1[sex]\n",
    "    #Compute correlation matrix\n",
    "    tempL = [item for sublist in [whtrDF.columns.tolist(), ['log_Base'+bwhtr]] for item in sublist]\n",
    "    tempDF1 = tempDF.drop(columns=tempL)\n",
    "    tempDF1 = tempDF1.corr(method='pearson')\n",
    "    print('• Combinations:', int(len(tempDF1)*(len(tempDF1)-1)/2))\n",
    "    \n",
    "    #Extract lower triangle matrix\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",
    "    print('• nrows:', len(tempDF1))\n",
    "    print('• |Pearson\\'s r| > 0.8:', len(tempDF1.loc[abs(tempDF1['Pearson_r'])>0.8]))\n",
    "    display(tempDF1.loc[abs(tempDF1['Pearson_r'])>0.8])\n",
    "    \n",
    "    tempD2[sex] = tempDF1\n",
    "\n",
    "#Distribution\n",
    "tempD3 = {'Female':'tab:red', 'Male':'tab:blue'}\n",
    "sns.set(style='ticks', font='Arial', context='notebook')\n",
    "plt.figure(figsize=(4, 3))\n",
    "for sex in tempD2.keys():\n",
    "    tempS = tempD2[sex]['Pearson_r']\n",
    "    sns.distplot(tempS, color=tempD3[sex], label=sex).set(xlim=(-1, 1), xticks=np.arange(-1, 1.1, 0.5))\n",
    "sns.despine()\n",
    "plt.xlabel('Pearson\\'s '+r'$r$')\n",
    "plt.ylabel('Density')\n",
    "plt.legend(loc='upper right')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Both sex cohort\n",
    "tempDF = chemDF_B\n",
    "bwhtr = 'ChemWHtR'\n",
    "bwhtr_color = 'g'\n",
    "\n",
    "print('Both sex cohort:')\n",
    "#Compute correlation matrix\n",
    "tempL = [item for sublist in [whtrDF.columns.tolist(), ['log_Base'+bwhtr]] for item in sublist]\n",
    "tempDF = tempDF.drop(columns=tempL)\n",
    "tempDF = tempDF.corr(method='pearson')\n",
    "print('• Combinations:', int(len(tempDF)*(len(tempDF)-1)/2))\n",
    "\n",
    "#Extract lower triangle matrix\n",
    "tempDF = tempDF.where(np.tril(np.ones(tempDF.shape), k=-1).astype(np.bool), other=np.nan)\n",
    "tempDF.index.rename('Variable1', inplace=True)\n",
    "tempDF = tempDF.reset_index().melt(var_name='Variable2', value_name='Pearson_r', id_vars=['Variable1'])\n",
    "tempDF = tempDF.dropna()\n",
    "print('• nrows:', len(tempDF))\n",
    "print('• |Pearson\\'s r| > 0.8:', len(tempDF.loc[abs(tempDF['Pearson_r'])>0.8]))\n",
    "display(tempDF.loc[abs(tempDF['Pearson_r'])>0.8])\n",
    "\n",
    "#Distribution\n",
    "sns.set(style='ticks', font='Arial', context='talk')\n",
    "plt.figure(figsize=(4, 3))\n",
    "sns.distplot(tempDF['Pearson_r'], color=bwhtr_color).set(xlim=(-1, 1), xticks=np.arange(-1, 1.1, 0.5))\n",
    "plt.axvline(x=tempDF['Pearson_r'].max(), **{'linestyle':'--', 'color':bwhtr_color})\n",
    "plt.axvline(x=tempDF['Pearson_r'].min(), **{'linestyle':'--', 'color':bwhtr_color})\n",
    "sns.despine()\n",
    "plt.ylabel('Density')\n",
    "plt.xlabel('Pearson\\'s '+r'$r$')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4-3. Prepare the beta-coefficients of LASSO models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bwhtr = 'ChemWHtR'\n",
    "tempD = {}\n",
    "for sex in ['Female', 'Male', 'BothSex']:\n",
    "    #Import the LASSO beta-coefficients\n",
    "    fileDir = './ExportData/'\n",
    "    ipynbName = '220822_Multiomics-BMI-NatMed1stRevision_WHtR-baseline-LASSO-ver2_'\n",
    "    fileName = bwhtr+'-'+sex+'-LASSObcoefs.tsv'\n",
    "    tempDF = pd.read_csv(fileDir+ipynbName+fileName, sep='\\t').set_index('Variable')\n",
    "    tempDF = tempDF.drop(index=['Intercept'])\n",
    "    tempD[sex] = tempDF\n",
    "    \n",
    "    #Check\n",
    "    print(sex+':')\n",
    "    print(' - Variables:', len(tempDF))\n",
    "    tempDF1 = tempDF.loc[tempDF['nZeros']!=10]\n",
    "    print(' - Variables with non-zero beta-coefficient:', len(tempDF1),\n",
    "          '(', len(tempDF1)/len(tempDF)*100, '%)')\n",
    "    tempDF1 = tempDF.loc[tempDF['nZeros']==0]\n",
    "    print(' - Variables with non-zero beta-coefficient in all 10 models:', len(tempDF1),\n",
    "          '(', len(tempDF1)/len(tempDF)*100, '%)')\n",
    "    tempDF1 = tempDF1.sort_values(by='Mean', ascending=False)\n",
    "    display(tempDF1)\n",
    "    print('')\n",
    "\n",
    "chemWHtR_F_bcoefs = tempD['Female']\n",
    "chemWHtR_M_bcoefs = tempD['Male']\n",
    "chemWHtR_B_bcoefs = tempD['BothSex']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4-4. Variables retained in all 10 models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempDF1 = chemWHtR_F_bcoefs\n",
    "tempDF2 = chemWHtR_M_bcoefs\n",
    "tempDF3 = chemWHtR_B_bcoefs\n",
    "category = 'Clinical labs'\n",
    "\n",
    "#Variables with non-zero beta-coefficient in all 10 models\n",
    "tempS1 = set(tempDF1.loc[tempDF1['nZeros']==0].index)\n",
    "tempS2 = set(tempDF2.loc[tempDF2['nZeros']==0].index)\n",
    "tempS3 = set(tempDF3.loc[tempDF3['nZeros']==0].index)\n",
    "\n",
    "#Common the robust variables\n",
    "var_111 = tempS1 & tempS2 & tempS3\n",
    "print('Common variables with non-zero beta-coefficient in all 10 models:\\n', var_111)\n",
    "var_110 = (tempS1 & tempS2) - var_111\n",
    "var_101 = (tempS1 & tempS3) - var_111\n",
    "var_011 = (tempS2 & tempS3) - var_111\n",
    "var_100 = tempS1 - (var_110 | var_101 | var_111)\n",
    "var_010 = tempS2 - (var_110 | var_011 | var_111)\n",
    "var_001 = tempS3 - (var_101 | var_011 | var_111)\n",
    "\n",
    "#The number of each subset\n",
    "var_subset = [len(var_100), len(var_010), len(var_110), len(var_001),\n",
    "              len(var_101), len(var_011), len(var_111)]\n",
    "print('Each subset:', var_subset)\n",
    "\n",
    "#Venn diagram\n",
    "sns.set(font='Arial', context='talk')\n",
    "venn3(subsets=var_subset, set_labels=('Female', 'Male', 'Both sex'),\n",
    "      set_colors=('tab:red', 'tab:blue', 'tab:green'), alpha=0.4)\n",
    "venn3_circles(subsets=var_subset)\n",
    "plt.title('Robust variables\\n—'+category+'—', fontdict={'fontsize':24})\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempDF1 = chemWHtR_F_bcoefs\n",
    "tempDF2 = chemWHtR_M_bcoefs\n",
    "tempDF3 = chemWHtR_B_bcoefs\n",
    "bwhtr_color = 'g'\n",
    "\n",
    "#Variables with non-zero beta-coefficient in all 10 models\n",
    "tempS1 = set(tempDF1.loc[tempDF1['nZeros']==0].index)\n",
    "tempS2 = set(tempDF2.loc[tempDF2['nZeros']==0].index)\n",
    "tempS3 = set(tempDF3.loc[tempDF3['nZeros']==0].index)\n",
    "\n",
    "#Union of the robust variables\n",
    "tempL = list(tempS1 | tempS2 | tempS3)\n",
    "print('Union of the variables with non-zero beta-coefficient in all 10 models:', len(tempL))\n",
    "\n",
    "#Merge (use original bcoefs to rescue ones with zeros in any 10 models)\n",
    "tempDF1 = tempDF1.loc[tempDF1.index.isin(tempL)].drop(columns=['Mean', 'SD', 'nZeros'])\n",
    "tempDF2 = tempDF2.loc[tempDF2.index.isin(tempL)].drop(columns=['Mean', 'SD', 'nZeros'])\n",
    "tempDF3 = tempDF3.loc[tempDF3.index.isin(tempL)].drop(columns=['Mean', 'SD', 'nZeros'])\n",
    "tempDF = pd.concat([tempDF1, tempDF2, tempDF3], axis=0)\n",
    "tempDF['Model'] = np.repeat(['Female', 'Male', 'Both sex'], [len(tempDF1), len(tempDF2), len(tempDF3)])\n",
    "\n",
    "#tidyr::gather in R\n",
    "tempDF1 = tempDF.reset_index().melt(var_name='LASSO', value_name='bcoef', id_vars=['Variable', 'Model'])\n",
    "\n",
    "#Order for plot\n",
    "tempDF2 = tempDF1.loc[tempDF1['Model']=='Both sex'].groupby('Variable', as_index=True).agg({'bcoef':np.mean})\n",
    "tempDF2 = tempDF2.sort_values(by='bcoef', ascending=False)\n",
    "\n",
    "#Plot\n",
    "sns.set(style='ticks', font='Arial', context='talk')\n",
    "sns.catplot(data=tempDF1, y='Variable', order=tempDF2.index.tolist(), x='bcoef', kind='box',\n",
    "            hue='Model', palette='Set1', dodge=True, height=9, aspect=1.1,\n",
    "            showfliers=True, flierprops={'marker':'o', 'markerfacecolor':'gray', 'alpha':0.4})\n",
    "plt.grid(axis='x', linestyle='--', color='black')\n",
    "for row_i in range(len(tempDF2)):\n",
    "    if row_i%2 == 0:\n",
    "        plt.axhspan(ymin=row_i-0.5, ymax=row_i+0.5, facecolor=bwhtr_color, alpha=0.2, zorder=0)\n",
    "plt.xlabel(r'$\\beta$'+'-coefficient in LASSO model')\n",
    "plt.ylabel('')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Both sex model\n",
    "tempDF = chemWHtR_B_bcoefs\n",
    "bwhtr_color = 'g'\n",
    "bwhtr = 'ChemWHtR'\n",
    "sex = 'BothSex'\n",
    "\n",
    "#Variables with non-zero beta-coefficient in all 10 models\n",
    "tempDF = tempDF.loc[tempDF['nZeros']==0].sort_values(by='Mean', ascending=False)\n",
    "tempDF = tempDF.drop(columns=['Mean', 'SD', 'nZeros'])\n",
    "print('Variables with non-zero beta-coefficient in all 10 models:', len(tempDF))\n",
    "\n",
    "#tidyr::gather in R\n",
    "tempDF1 = tempDF.reset_index().melt(var_name='LASSO', value_name='bcoef', id_vars=['Variable'])\n",
    "\n",
    "#Plot\n",
    "sns.set(style='ticks', font='Arial', context='talk')\n",
    "plt.figure(figsize=(5, 8))\n",
    "p = sns.boxplot(data=tempDF1, y='Variable', x='bcoef', color=bwhtr_color, dodge=False, saturation=1,\n",
    "                showfliers=True, flierprops={'marker':'o', 'markerfacecolor':'gray', 'alpha':0.4},\n",
    "                showcaps=True, notch=False)\n",
    "p.set(xlim=(-0.046, 0.046), xticks=np.arange(-0.04, 0.041, 0.02))#Fixed across omics\n",
    "p.grid(axis='x', linestyle='--', color='black')\n",
    "sns.despine()\n",
    "##Change default dull line color of sns.boxplot (saturation parameter is for patch)\n",
    "for line in p.get_lines():\n",
    "    line.set_color('k')\n",
    "for box in p.artists:\n",
    "    box.set_edgecolor('k')\n",
    "##Add background color\n",
    "for row_i in range(len(tempDF)):\n",
    "    if row_i%2 == 0:\n",
    "        plt.axhspan(ymin=row_i-0.5, ymax=row_i+0.5, facecolor=bwhtr_color, alpha=0.2, zorder=0)\n",
    "plt.ylabel('')\n",
    "plt.xlabel(r'$\\beta$'+'-coefficient in LASSO model\\n [log-scaled WHtR per s.d.]')\n",
    "##Save\n",
    "fileDir = './ExportFigures/'\n",
    "ipynbName = '220823_Multiomics-BMI-NatMed1stRevision_WHtR-LASSO-bcoef-ver2_'\n",
    "fileName = bwhtr+'-'+sex+'-bcoef_non-zero-in-all.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": [
    "### 4-5. Correlation b/w pairwise variables retained in all 10 models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Both sex model\n",
    "tempDF1 = chemWHtR_B_bcoefs\n",
    "tempDF2 = chemDF_B\n",
    "bwhtr_color = 'g'\n",
    "\n",
    "print('Both sex model:')\n",
    "#Variables with non-zero beta-coefficient in all 10 models\n",
    "tempL = tempDF1.loc[tempDF1['nZeros']==0].index.tolist()\n",
    "\n",
    "#Compute correlation matrix\n",
    "tempDF = tempDF2[tempL]\n",
    "tempDF = tempDF.corr(method='pearson')\n",
    "print('• Combinations:', int(len(tempDF)*(len(tempDF)-1)/2))\n",
    "\n",
    "#Extract lower triangle matrix\n",
    "tempDF = tempDF.where(np.tril(np.ones(tempDF.shape), k=-1).astype(np.bool), other=np.nan)\n",
    "tempDF.index.rename('Variable1', inplace=True)\n",
    "tempDF = tempDF.reset_index().melt(var_name='Variable2', value_name='Pearson_r', id_vars=['Variable1'])\n",
    "tempDF = tempDF.dropna()\n",
    "print('• nrows:', len(tempDF))\n",
    "print('• |Pearson\\'s r| > 0.8:', len(tempDF.loc[abs(tempDF['Pearson_r'])>0.8]))\n",
    "display(tempDF.loc[abs(tempDF['Pearson_r'])>0.8])\n",
    "\n",
    "#Distribution\n",
    "sns.set(style='ticks', font='Arial', context='talk')\n",
    "plt.figure(figsize=(4, 3))\n",
    "sns.distplot(tempDF['Pearson_r'], color=bwhtr_color).set(xlim=(-1, 1), xticks=np.arange(-1, 1.1, 0.5))\n",
    "plt.axvline(x=tempDF['Pearson_r'].max(), **{'linestyle':'--', 'color':bwhtr_color})\n",
    "plt.axvline(x=tempDF['Pearson_r'].min(), **{'linestyle':'--', 'color':bwhtr_color})\n",
    "sns.despine()\n",
    "plt.ylabel('Density')\n",
    "plt.xlabel('Pearson\\'s '+r'$r$')\n",
    "plt.show()\n",
    "\n",
    "print('Variables with non-zero beta-coefficient in all 10 models:', len(tempL))\n",
    "#Clustermap\n",
    "tempDF = tempDF2[tempL]\n",
    "tempDF = tempDF.corr(method='pearson')\n",
    "sns.set(style='ticks', font='Arial', context='notebook')\n",
    "cmap = sns.diverging_palette(220, 20, as_cmap=True)\n",
    "cm = sns.clustermap(tempDF, method='ward', metric='euclidean', cmap=cmap,\n",
    "                    row_cluster=True, col_cluster=True, row_linkage=None, col_linkage=None,\n",
    "                    row_colors=None, col_colors=None, xticklabels=True, yticklabels=True,\n",
    "                    dendrogram_ratio=(0.15, 0.15), cbar_pos=(0.85, 0.07, 0.02, 0.1),\n",
    "                    figsize=(7, 7), **{'center':0, 'vmin':-1, 'vmax':1})\n",
    "cm.cax.set_title('Pearson\\'s '+r'$r$')\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.5])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4-6. Explained variance in WHtR by the variables retained in ≥1 model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> OLS linear regression model for significance: log(WHtR) ~ b0 + b1\\*Variable + b2\\*Sex + b3\\*BaseAge + b4\\*PCs  \n",
    "> Univariate model: log(WHtR) ~ b0 + b1\\*Variable  \n",
    ">\n",
    "> Of note, biological WHtR is included for comparison.  \n",
    "> As well as LASSO models, unstandardized WHtR is used for the dependent variable.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempD1 = {'Female':chemWHtR_F_bcoefs, 'Male':chemWHtR_M_bcoefs, 'BothSex':chemWHtR_B_bcoefs}\n",
    "tempD2 = {'Female':chemDF_F, 'Male':chemDF_M, 'BothSex':chemDF_B}\n",
    "tempD3 = {'Female':whtrDF_F, 'Male':whtrDF_M, 'BothSex':whtrDF_B}\n",
    "yvar = 'WHtR'\n",
    "bwhtr = 'ChemWHtR'\n",
    "topX = 30\n",
    "for sex in tempD1.keys():\n",
    "    print(sex+' model:')\n",
    "    \n",
    "    #Variables with non-zero beta-coefficient at least 1 of 10 models\n",
    "    tempDF = tempD1[sex]\n",
    "    tempL = tempDF.loc[tempDF['nZeros']==10].index.tolist()#To be excluded\n",
    "    \n",
    "    #Prepare DF for regressions\n",
    "    tempDF = tempD2[sex].drop(columns=tempL)\n",
    "    tempDF = tempDF.rename(columns={'log_Base'+bwhtr:bwhtr+' model'})\n",
    "    #Replace the standardized dependent variable with the unstandardized\n",
    "    tempDF1 = tempD3[sex]\n",
    "    tempDF['log_Base'+yvar] = tempDF1['log_Base'+yvar]\n",
    "    \n",
    "    #Independent variables for regressions (including bWHtR)\n",
    "    tempL = tempDF.drop(columns=whtrDF.columns.tolist()).columns.tolist()\n",
    "    \n",
    "    #Perform OLS linear regression\n",
    "    tempL1 = []#For sample size\n",
    "    tempL2 = []#For R2 in the univariate model\n",
    "    tempL3 = []#For R2\n",
    "    tempL4 = []#For beta-coefficient\n",
    "    tempL5 = []#For SE of beta-coefficient\n",
    "    tempL6 = []#For t-statistic\n",
    "    tempL7 = []#For residual degrees of freedom\n",
    "    tempL8 = []#For P-value\n",
    "    t_start = time.time()\n",
    "    for var in tempL:\n",
    "        #Rename independent variable\n",
    "        tempDF1 = tempDF.rename(columns={var: 'Variable'})\n",
    "        #Add a constant for the intercept -> Similar to R, smf automatically add a constant\n",
    "        \n",
    "        #Save sample size\n",
    "        tempL1.append(len(tempDF1))\n",
    "        \n",
    "        #Fit univariate model\n",
    "        formula = 'log_Base'+yvar+' ~ Variable'\n",
    "        fit_res = smf.ols(formula, data=tempDF1).fit()\n",
    "        #Save R2 [%]\n",
    "        tempL2.append(fit_res.rsquared*100)\n",
    "\n",
    "        #Fit full model\n",
    "        if sex=='BothSex':\n",
    "            formula = 'log_Base'+yvar+' ~ Variable + C(Sex) + BaseAge + PC1 + PC2 + PC3 + PC4 + PC5'\n",
    "        else:#Sex is invariant\n",
    "            formula = 'log_Base'+yvar+' ~ Variable + BaseAge + PC1 + PC2 + PC3 + PC4 + PC5'\n",
    "        fit_res = smf.ols(formula, data=tempDF1).fit()\n",
    "        #Save R2 [%]\n",
    "        tempL3.append(fit_res.rsquared*100)\n",
    "        #Save beta-coefficient of the variable\n",
    "        tempL4.append(fit_res.params['Variable'])\n",
    "        tempL5.append(fit_res.bse['Variable'])\n",
    "        #Save t-statistic of the variable\n",
    "        tempL6.append(fit_res.tvalues['Variable'])\n",
    "        #Save residual degrees of freedom\n",
    "        tempL7.append(int(fit_res.df_resid))\n",
    "        #Save P-value of the variable\n",
    "        tempL8.append(fit_res.pvalues['Variable'])\n",
    "    t_elapsed = time.time() - t_start\n",
    "    print('Elapsed time for', len(tempL), 'OLS linear regressions (including '+bwhtr+'):',\n",
    "          round(t_elapsed//60), 'min', round(t_elapsed%60, 1), 'sec')\n",
    "    \n",
    "    #Clean the results\n",
    "    tempDF = pd.DataFrame({'N':tempL1, 'UnivarR2':tempL2, 'R2':tempL3,\n",
    "                           'OLSbcoef':tempL4, 'OLSbcoefSE':tempL5,\n",
    "                           'tStat':tempL6, 'DoF':tempL7,  'Pval':tempL8},\n",
    "                          index=pd.Index(tempL, name='Variable'))\n",
    "    ##P-value adjustment by using Benjamini–Hochberg method\n",
    "    tempDF['AdjPval'] = multi.multipletests(tempDF['Pval'], alpha=0.05, method='fdr_bh',\n",
    "                                            is_sorted=False, returnsorted=False)[1]\n",
    "    ##Add the LASSO results\n",
    "    tempDF1 = tempD1[sex]\n",
    "    tempDF['LASSObcoef'] = tempDF1.loc[tempDF1.index.isin(tempL)]['Mean']#NaN in bWHtR reference\n",
    "    tempDF['LASSOnZeros'] = tempDF1.loc[tempDF1.index.isin(tempL)]['nZeros']#NaN in bWHtR reference\n",
    "    tempDF = tempDF.sort_values(by='AdjPval', ascending=True)\n",
    "    \n",
    "    #Save the cleaned DF\n",
    "    fileDir = './ExportData/'\n",
    "    ipynbName = '220823_Multiomics-BMI-NatMed1stRevision_WHtR-LASSO-bcoef-ver2_'\n",
    "    fileName = bwhtr+'-'+sex+'-OLS.tsv'\n",
    "    tempDF.to_csv(fileDir+ipynbName+fileName, sep='\\t', index=True)\n",
    "    \n",
    "    #Extact significant variables\n",
    "    tempDF = tempDF[tempDF['AdjPval']<0.05]\n",
    "    print('Variables significantly associated with WHtR (FDR < 0.05):', len(tempDF)-1)#bWHtR reference\n",
    "    ##Top X (+ bWHtR)\n",
    "    tempDF = tempDF.iloc[0:(topX+1)]\n",
    "    display(tempDF)\n",
    "    ##Category and color\n",
    "    tempL = []\n",
    "    for row_n in tempDF.index.tolist():\n",
    "        if row_n==bwhtr+' model':\n",
    "            tempL.append('Positive association')\n",
    "        else:\n",
    "            bcoef = tempDF.loc[row_n, 'OLSbcoef']\n",
    "            if bcoef>0:\n",
    "                tempL.append('Positive association')\n",
    "            elif bcoef<0:\n",
    "                tempL.append('Negative association')\n",
    "            else:\n",
    "                tempL.append('Error?')\n",
    "    tempDF['Association'] = tempL\n",
    "    tempD = {'Positive association':'tab:red', 'Negative association':'tab:blue'}\n",
    "    ##Plot\n",
    "    tempDF = tempDF.reset_index()\n",
    "    sns.set(style='ticks', font='Arial', context='talk')\n",
    "    plt.figure(figsize=(10, 4))\n",
    "    p = sns.barplot(data=tempDF, y='UnivarR2', x='Variable', hue='Association', hue_order=tempD.keys(),\n",
    "                    dodge=False, palette=tempD, edgecolor='k')\n",
    "    p.set(ylim=(0, 80), yticks=np.arange(0, 80, 15))#Fixed across omics\n",
    "    p.grid(axis='y', linestyle='--', color='black')\n",
    "    sns.despine()\n",
    "    plt.xticks(rotation=90, horizontalalignment='center')\n",
    "    plt.xlabel('')\n",
    "    plt.ylabel('Explained variance in WHtR [%]')\n",
    "    plt.legend(loc='upper right')\n",
    "    if sex=='BothSex':#Save\n",
    "        fileDir = './ExportFigures/'\n",
    "        ipynbName = '220823_Multiomics-BMI-NatMed1stRevision_WHtR-LASSO-bcoef-ver2_'\n",
    "        fileName = bwhtr+'-'+sex+'-OLS-top'+str(topX)+'.tif'\n",
    "        plt.gcf().savefig(fileDir+ipynbName+fileName, dpi=300, bbox_inches='tight', pad_inches=0.04,\n",
    "                          pil_kwargs={'compression':'tiff_lzw'})\n",
    "    plt.show()\n",
    "    print('')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Metabolomics, Proteomics, and Clinical labs-combined omics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5-1. Prepare and standardize variables"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Variables, including bWHtR (as a reference) and covariates, are standardized for OLS linear regression.  \n",
    "> –> In this version, unstandardized WHtR is used as the dependent variable. However, WHtR within the following DFs is standardized (just for the code simplicity).  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempDF = combiDF\n",
    "tempD1 = {'Female':whtrDF_F, 'Male':whtrDF_M, 'BothSex':whtrDF_B}\n",
    "bwhtr = 'CombiWHtR'\n",
    "tempD2 = {}\n",
    "for sex in tempD1.keys():\n",
    "    #Add covariates (and WHtR) while selecting the cohort\n",
    "    tempDF1 = tempD1[sex]\n",
    "    tempDF1 = pd.merge(tempDF, tempDF1, left_index=True, right_index=True, how='inner')\n",
    "    \n",
    "    #Add biological WHtR\n",
    "    fileDir = './ExportData/'\n",
    "    ipynbName = '220822_Multiomics-BMI-NatMed1stRevision_WHtR-baseline-LASSO-ver2_'\n",
    "    fileName = bwhtr+'-'+sex+'.tsv'\n",
    "    tempDF2 = pd.read_csv(fileDir+ipynbName+fileName, sep='\\t', dtype={'public_client_id':str})\n",
    "    tempDF2 = tempDF2.set_index('public_client_id')\n",
    "    tempDF1['log_Base'+bwhtr] = tempDF2['log_Base'+bwhtr]\n",
    "    \n",
    "    #Z-score transformation\n",
    "    tempDF2 = tempDF1.select_dtypes(include=[np.number])\n",
    "    scaler = StandardScaler(copy=True, with_mean=True, with_std=True)\n",
    "    tempA = scaler.fit_transform(tempDF2)#Column direction\n",
    "    tempDF2 = pd.DataFrame(data=tempA, index=tempDF2.index, columns=tempDF2.columns)\n",
    "    ##Recover categorical covariates\n",
    "    tempDF3 = tempDF1.select_dtypes(exclude=[np.number])\n",
    "    tempDF1 = pd.merge(tempDF2, tempDF3, left_index=True, right_index=True, how='left')\n",
    "    \n",
    "    print(sex)\n",
    "    display(tempDF1.describe(include='all'))\n",
    "    tempD2[sex] = tempDF1\n",
    "\n",
    "combiDF_F = tempD2['Female']\n",
    "combiDF_M = tempD2['Male']\n",
    "combiDF_B = tempD2['BothSex']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5-2. Correlation b/w all pairwise variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Female and male cohorts\n",
    "tempD1 = {'Female':combiDF_F, 'Male':combiDF_M}\n",
    "bwhtr = 'CombiWHtR'\n",
    "tempD2 = {}\n",
    "for sex in tempD1.keys():\n",
    "    print(sex+' cohort:')\n",
    "    tempDF = tempD1[sex]\n",
    "    #Compute correlation matrix\n",
    "    tempL = [item for sublist in [whtrDF.columns.tolist(), ['log_Base'+bwhtr]] for item in sublist]\n",
    "    tempDF1 = tempDF.drop(columns=tempL)\n",
    "    tempDF1 = tempDF1.corr(method='pearson')\n",
    "    print('• Combinations:', int(len(tempDF1)*(len(tempDF1)-1)/2))\n",
    "    \n",
    "    #Extract lower triangle matrix\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",
    "    print('• nrows:', len(tempDF1))\n",
    "    print('• |Pearson\\'s r| > 0.8:', len(tempDF1.loc[abs(tempDF1['Pearson_r'])>0.8]))\n",
    "    display(tempDF1.loc[abs(tempDF1['Pearson_r'])>0.8])\n",
    "    \n",
    "    tempD2[sex] = tempDF1\n",
    "\n",
    "#Distribution\n",
    "tempD3 = {'Female':'tab:red', 'Male':'tab:blue'}\n",
    "sns.set(style='ticks', font='Arial', context='notebook')\n",
    "plt.figure(figsize=(4, 3))\n",
    "for sex in tempD2.keys():\n",
    "    tempS = tempD2[sex]['Pearson_r']\n",
    "    sns.distplot(tempS, color=tempD3[sex], label=sex).set(xlim=(-1, 1), xticks=np.arange(-1, 1.1, 0.5))\n",
    "sns.despine()\n",
    "plt.xlabel('Pearson\\'s '+r'$r$')\n",
    "plt.ylabel('Density')\n",
    "plt.legend(loc='upper right')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Both sex cohort\n",
    "tempDF = combiDF_B\n",
    "bwhtr = 'CombiWHtR'\n",
    "bwhtr_color = 'm'\n",
    "\n",
    "print('Both sex cohort:')\n",
    "#Compute correlation matrix\n",
    "tempL = [item for sublist in [whtrDF.columns.tolist(), ['log_Base'+bwhtr]] for item in sublist]\n",
    "tempDF = tempDF.drop(columns=tempL)\n",
    "tempDF = tempDF.corr(method='pearson')\n",
    "print('• Combinations:', int(len(tempDF)*(len(tempDF)-1)/2))\n",
    "\n",
    "#Extract lower triangle matrix\n",
    "tempDF = tempDF.where(np.tril(np.ones(tempDF.shape), k=-1).astype(np.bool), other=np.nan)\n",
    "tempDF.index.rename('Variable1', inplace=True)\n",
    "tempDF = tempDF.reset_index().melt(var_name='Variable2', value_name='Pearson_r', id_vars=['Variable1'])\n",
    "tempDF = tempDF.dropna()\n",
    "print('• nrows:', len(tempDF))\n",
    "print('• |Pearson\\'s r| > 0.8:', len(tempDF.loc[abs(tempDF['Pearson_r'])>0.8]))\n",
    "display(tempDF.loc[abs(tempDF['Pearson_r'])>0.8])\n",
    "\n",
    "#Distribution\n",
    "sns.set(style='ticks', font='Arial', context='talk')\n",
    "plt.figure(figsize=(4, 3))\n",
    "sns.distplot(tempDF['Pearson_r'], color=bwhtr_color).set(xlim=(-1, 1), xticks=np.arange(-1, 1.1, 0.5))\n",
    "plt.axvline(x=tempDF['Pearson_r'].max(), **{'linestyle':'--', 'color':bwhtr_color})\n",
    "plt.axvline(x=tempDF['Pearson_r'].min(), **{'linestyle':'--', 'color':bwhtr_color})\n",
    "sns.despine()\n",
    "plt.ylabel('Density')\n",
    "plt.xlabel('Pearson\\'s '+r'$r$')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5-3. Prepare the beta-coefficients of LASSO models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bwhtr = 'CombiWHtR'\n",
    "tempD = {}\n",
    "for sex in ['Female', 'Male', 'BothSex']:\n",
    "    #Import the LASSO beta-coefficients\n",
    "    fileDir = './ExportData/'\n",
    "    ipynbName = '220822_Multiomics-BMI-NatMed1stRevision_WHtR-baseline-LASSO-ver2_'\n",
    "    fileName = bwhtr+'-'+sex+'-LASSObcoefs.tsv'\n",
    "    tempDF = pd.read_csv(fileDir+ipynbName+fileName, sep='\\t').set_index('Variable')\n",
    "    tempDF = tempDF.drop(index=['Intercept'])\n",
    "    tempD[sex] = tempDF\n",
    "    \n",
    "    #Check\n",
    "    print(sex+':')\n",
    "    print(' - Variables:', len(tempDF))\n",
    "    tempDF1 = tempDF.loc[tempDF['nZeros']!=10]\n",
    "    print(' - Variables with non-zero beta-coefficient:', len(tempDF1),\n",
    "          '(', len(tempDF1)/len(tempDF)*100, '%)')\n",
    "    tempDF1 = tempDF.loc[tempDF['nZeros']==0]\n",
    "    print(' - Variables with non-zero beta-coefficient in all 10 models:', len(tempDF1),\n",
    "          '(', len(tempDF1)/len(tempDF)*100, '%)')\n",
    "    tempDF1 = tempDF1.sort_values(by='Mean', ascending=False)\n",
    "    display(tempDF1)\n",
    "    print('')\n",
    "\n",
    "combiWHtR_F_bcoefs = tempD['Female']\n",
    "combiWHtR_M_bcoefs = tempD['Male']\n",
    "combiWHtR_B_bcoefs = tempD['BothSex']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5-4. Variables retained in all 10 models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempDF1 = combiWHtR_F_bcoefs\n",
    "tempDF2 = combiWHtR_M_bcoefs\n",
    "tempDF3 = combiWHtR_B_bcoefs\n",
    "category = 'Combined omics'\n",
    "\n",
    "#Variables with non-zero beta-coefficient in all 10 models\n",
    "tempS1 = set(tempDF1.loc[tempDF1['nZeros']==0].index)\n",
    "tempS2 = set(tempDF2.loc[tempDF2['nZeros']==0].index)\n",
    "tempS3 = set(tempDF3.loc[tempDF3['nZeros']==0].index)\n",
    "\n",
    "#Common the robust variables\n",
    "var_111 = tempS1 & tempS2 & tempS3\n",
    "print('Common variables with non-zero beta-coefficient in all 10 models:\\n', var_111)\n",
    "var_110 = (tempS1 & tempS2) - var_111\n",
    "var_101 = (tempS1 & tempS3) - var_111\n",
    "var_011 = (tempS2 & tempS3) - var_111\n",
    "var_100 = tempS1 - (var_110 | var_101 | var_111)\n",
    "var_010 = tempS2 - (var_110 | var_011 | var_111)\n",
    "var_001 = tempS3 - (var_101 | var_011 | var_111)\n",
    "\n",
    "#The number of each subset\n",
    "var_subset = [len(var_100), len(var_010), len(var_110), len(var_001),\n",
    "              len(var_101), len(var_011), len(var_111)]\n",
    "print('Each subset:', var_subset)\n",
    "\n",
    "#Venn diagram\n",
    "sns.set(font='Arial', context='talk')\n",
    "venn3(subsets=var_subset, set_labels=('Female', 'Male', 'Both sex'),\n",
    "      set_colors=('tab:red', 'tab:blue', 'tab:green'), alpha=0.4)\n",
    "venn3_circles(subsets=var_subset)\n",
    "plt.title('Robust variables\\n—'+category+'—', fontdict={'fontsize':24})\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempDF1 = combiWHtR_F_bcoefs\n",
    "tempDF2 = combiWHtR_M_bcoefs\n",
    "tempDF3 = combiWHtR_B_bcoefs\n",
    "\n",
    "#Variables with non-zero beta-coefficient in all 10 models\n",
    "tempS1 = set(tempDF1.loc[tempDF1['nZeros']==0].index)\n",
    "tempS2 = set(tempDF2.loc[tempDF2['nZeros']==0].index)\n",
    "tempS3 = set(tempDF3.loc[tempDF3['nZeros']==0].index)\n",
    "\n",
    "#Union of the robust variables\n",
    "tempL = list(tempS1 | tempS2 | tempS3)\n",
    "print('Union of the variables with non-zero beta-coefficient in all 10 models:', len(tempL))\n",
    "\n",
    "#Merge (use original bcoefs to rescue ones with zeros in any 10 models)\n",
    "tempDF1 = tempDF1.loc[tempDF1.index.isin(tempL)].drop(columns=['Mean', 'SD', 'nZeros'])\n",
    "tempDF2 = tempDF2.loc[tempDF2.index.isin(tempL)].drop(columns=['Mean', 'SD', 'nZeros'])\n",
    "tempDF3 = tempDF3.loc[tempDF3.index.isin(tempL)].drop(columns=['Mean', 'SD', 'nZeros'])\n",
    "tempDF = pd.concat([tempDF1, tempDF2, tempDF3], axis=0)\n",
    "tempDF['Model'] = np.repeat(['Female', 'Male', 'Both sex'], [len(tempDF1), len(tempDF2), len(tempDF3)])\n",
    "\n",
    "#tidyr::gather in R\n",
    "tempDF1 = tempDF.reset_index().melt(var_name='LASSO', value_name='bcoef', id_vars=['Variable', 'Model'])\n",
    "\n",
    "#Order for plot and category for shading\n",
    "tempDF2 = tempDF1.loc[tempDF1['Model']=='Both sex'].groupby('Variable', as_index=True).agg({'bcoef':np.mean})\n",
    "tempDF2 = tempDF2.sort_values(by='bcoef', ascending=False)\n",
    "tempDF2 = pd.DataFrame(index=tempDF2.index)\n",
    "tempL1 = []\n",
    "tempL2 = []\n",
    "for row_n in tempDF2.index.tolist():\n",
    "    if row_n in metDF.columns.tolist():\n",
    "        tempL1.append('Metabolomics')\n",
    "        tempL2.append('b')\n",
    "    elif row_n in chemDF.columns.tolist():\n",
    "        tempL1.append('Clinical labs')\n",
    "        tempL2.append('g')\n",
    "    elif row_n in protDF.columns.tolist():\n",
    "        tempL1.append('Proteomics')\n",
    "        tempL2.append('r')\n",
    "    else:\n",
    "        tempL1.append('Error?')\n",
    "        tempL2.append('k')\n",
    "tempDF2['Category'] = tempL1\n",
    "tempDF2['Color'] = tempL2\n",
    "\n",
    "#Plot\n",
    "sns.set(style='ticks', font='Arial', context='talk')\n",
    "sns.catplot(data=tempDF1, y='Variable', order=tempDF2.index.tolist(), x='bcoef', kind='box',\n",
    "            hue='Model', palette='Set1', dodge=True, height=20, aspect=0.7,\n",
    "            showfliers=True, flierprops={'marker':'o', 'markerfacecolor':'gray', 'alpha':0.4})\n",
    "plt.grid(axis='x', linestyle='--', color='black')\n",
    "for row_i in range(len(tempDF2)):\n",
    "    cat_color = tempDF2['Color'].iloc[row_i]\n",
    "    if row_i%2 == 0:\n",
    "        plt.axhspan(ymin=row_i-0.5, ymax=row_i+0.5, facecolor=cat_color, alpha=0.2, zorder=0)\n",
    "    else:\n",
    "        plt.axhspan(ymin=row_i-0.5, ymax=row_i+0.5, facecolor=cat_color, alpha=0.4, zorder=0)\n",
    "plt.xlabel(r'$\\beta$'+'-coefficient in LASSO model')\n",
    "plt.ylabel('')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Both sex model\n",
    "tempDF = combiWHtR_B_bcoefs\n",
    "bwhtr = 'CombiWHtR'\n",
    "sex = 'BothSex'\n",
    "\n",
    "#Variables with non-zero beta-coefficient in all 10 models\n",
    "tempDF = tempDF.loc[tempDF['nZeros']==0].sort_values(by='Mean', ascending=False)\n",
    "tempDF = tempDF.drop(columns=['Mean', 'SD', 'nZeros'])\n",
    "print('Variables with non-zero beta-coefficient in all 10 models:', len(tempDF))\n",
    "\n",
    "#tidyr::gather in R\n",
    "tempDF1 = tempDF.reset_index().melt(var_name='LASSO', value_name='bcoef', id_vars=['Variable'])\n",
    "\n",
    "#Category for shading\n",
    "tempDF2 = pd.DataFrame(index=tempDF.index)\n",
    "tempL1 = []\n",
    "tempL2 = []\n",
    "for row_n in tempDF2.index.tolist():\n",
    "    if row_n in metDF.columns.tolist():\n",
    "        tempL1.append('Metabolomics')\n",
    "        tempL2.append('b')\n",
    "    elif row_n in chemDF.columns.tolist():\n",
    "        tempL1.append('Clinical labs')\n",
    "        tempL2.append('g')\n",
    "    elif row_n in protDF.columns.tolist():\n",
    "        tempL1.append('Proteomics')\n",
    "        tempL2.append('r')\n",
    "    else:\n",
    "        tempL1.append('Error?')\n",
    "        tempL2.append('k')\n",
    "tempDF2['Category'] = tempL1\n",
    "tempDF2['Color'] = tempL2\n",
    "\n",
    "#Plot\n",
    "sns.set(style='ticks', font='Arial', context='talk')\n",
    "plt.figure(figsize=(5, 12.5))\n",
    "p = sns.boxplot(data=tempDF1, y='Variable', order=tempDF2.index.tolist(), x='bcoef',\n",
    "                palette=tempDF2['Color'], dodge=False, saturation=1,\n",
    "                showfliers=True, flierprops={'marker':'o', 'markerfacecolor':'gray', 'alpha':0.4},\n",
    "                showcaps=True, notch=False)\n",
    "#p.set(xlim=(-0.046, 0.046), xticks=np.arange(-0.04, 0.041, 0.02))#Fixed across omics\n",
    "p.grid(axis='x', linestyle='--', color='k')\n",
    "sns.despine()\n",
    "##Change default dull line color of sns.boxplot (saturation parameter is for patch)\n",
    "for line in p.get_lines():\n",
    "    line.set_color('k')\n",
    "for box in p.artists:\n",
    "    box.set_edgecolor('k')\n",
    "##Add background color\n",
    "for row_i in range(len(tempDF2)):\n",
    "    cat_color = tempDF2['Color'].iloc[row_i]\n",
    "    if row_i%2 == 0:\n",
    "        plt.axhspan(ymin=row_i-0.5, ymax=row_i+0.5, facecolor=cat_color, alpha=0.4, zorder=0)\n",
    "    else:\n",
    "        plt.axhspan(ymin=row_i-0.5, ymax=row_i+0.5, facecolor=cat_color, alpha=0.4, zorder=0)\n",
    "plt.ylabel('')\n",
    "plt.xlabel(r'$\\beta$'+'-coefficient in LASSO model\\n [log-scaled WHtR per s.d.]')\n",
    "#Add legend\n",
    "legend1 = mpatches.Patch(facecolor='b', edgecolor='k', label='Metabolomics')\n",
    "legend2 = mpatches.Patch(facecolor='r', edgecolor='k', label='Proteomics')\n",
    "legend3 = mpatches.Patch(facecolor='g', edgecolor='k', label='Clinical labs')\n",
    "plt.legend(handles=[legend1, legend2, legend3], fontsize='large',\n",
    "           title='Omics category', title_fontsize='x-large',\n",
    "           bbox_to_anchor=(-0.7, 0.0), loc='lower right', borderaxespad=0)#Manual adjustment\n",
    "##Save\n",
    "fileDir = './ExportFigures/'\n",
    "ipynbName = '220823_Multiomics-BMI-NatMed1stRevision_WHtR-LASSO-bcoef-ver2_'\n",
    "fileName = bwhtr+'-'+sex+'-bcoef_non-zero-in-all.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": [
    "### 5-5. Correlation b/w pairwise variables retained in all 10 models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Both sex model\n",
    "tempDF1 = combiWHtR_B_bcoefs\n",
    "tempDF2 = combiDF_B\n",
    "bwhtr_color = 'm'\n",
    "bwhtr = 'CombiWHtR'\n",
    "sex = 'BothSex'\n",
    "\n",
    "print('Both sex model:')\n",
    "#Variables with non-zero beta-coefficient in all 10 models\n",
    "tempL = tempDF1.loc[tempDF1['nZeros']==0].index.tolist()\n",
    "\n",
    "#Compute correlation matrix\n",
    "tempDF = tempDF2[tempL]\n",
    "tempDF = tempDF.corr(method='pearson')\n",
    "print('• Combinations:', int(len(tempDF)*(len(tempDF)-1)/2))\n",
    "\n",
    "#Extract lower triangle matrix\n",
    "tempDF = tempDF.where(np.tril(np.ones(tempDF.shape), k=-1).astype(np.bool), other=np.nan)\n",
    "tempDF.index.rename('Variable1', inplace=True)\n",
    "tempDF = tempDF.reset_index().melt(var_name='Variable2', value_name='Pearson_r', id_vars=['Variable1'])\n",
    "tempDF = tempDF.dropna()\n",
    "print('• nrows:', len(tempDF))\n",
    "print('• |Pearson\\'s r| > 0.8:', len(tempDF.loc[abs(tempDF['Pearson_r'])>0.8]))\n",
    "display(tempDF.loc[abs(tempDF['Pearson_r'])>0.8])\n",
    "\n",
    "#Distribution\n",
    "sns.set(style='ticks', font='Arial', context='talk')\n",
    "plt.figure(figsize=(4, 3))\n",
    "sns.distplot(tempDF['Pearson_r'], color=bwhtr_color).set(xlim=(-1, 1), xticks=np.arange(-1, 1.1, 0.5))\n",
    "plt.axvline(x=tempDF['Pearson_r'].max(), **{'linestyle':'--', 'color':bwhtr_color})\n",
    "plt.axvline(x=tempDF['Pearson_r'].min(), **{'linestyle':'--', 'color':bwhtr_color})\n",
    "sns.despine()\n",
    "plt.ylabel('Density')\n",
    "plt.xlabel('Pearson\\'s '+r'$r$')\n",
    "plt.show()\n",
    "\n",
    "print('Variables with non-zero beta-coefficient in all 10 models:', len(tempL))\n",
    "tempDF = tempDF2[tempL]\n",
    "tempDF = tempDF.corr(method='pearson')\n",
    "#Category\n",
    "tempS = pd.Series(index=tempDF.index, name='Omics category')\n",
    "for row_n in tempS.index.tolist():\n",
    "    if row_n in metDF.columns.tolist():\n",
    "        tempS[row_n] = 'tab:blue'\n",
    "    elif row_n in chemDF.columns.tolist():\n",
    "        tempS[row_n] = 'tab:green'\n",
    "    elif row_n in protDF.columns.tolist():\n",
    "        tempS[row_n] = 'tab:red'\n",
    "    else:#Just in case\n",
    "        tempS[row_n] = 'k'\n",
    "#Clustermap\n",
    "sns.set(style='ticks', font='Arial', context='notebook')\n",
    "cmap = sns.diverging_palette(220, 20, as_cmap=True)\n",
    "cm = sns.clustermap(tempDF, method='ward', metric='euclidean', cmap=cmap,\n",
    "                    row_cluster=True, col_cluster=True, row_linkage=None, col_linkage=None,\n",
    "                    row_colors=tempS, col_colors=tempS, xticklabels=True, yticklabels=True,\n",
    "                    dendrogram_ratio=(0.1, 0.1), colors_ratio=(0.02, 0.02),\n",
    "                    cbar_pos=(0.7, 0.2, 0.25, 0.025), cbar_kws={\"orientation\": \"horizontal\"},\n",
    "                    figsize=(11, 11), **{'center':0, 'vmin':-1, 'vmax':1})\n",
    "cm.cax.set_title('Pearson\\'s '+r'$r$', size='xx-large', verticalalignment='bottom')\n",
    "cm.cax.tick_params(labelsize='x-large')\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.5])\n",
    "##row/column color bar legend (axis is same with cm.cax!)\n",
    "legend1 = mpatches.Patch(color='tab:blue', label='Metabolomics')\n",
    "legend2 = mpatches.Patch(color='tab:red', label='Proteomics')\n",
    "legend3 = mpatches.Patch(color='tab:green', label='Clinical labs')\n",
    "plt.legend(handles=[legend1, legend2, legend3], fontsize='x-large',\n",
    "           title=tempS.name, title_fontsize='xx-large',\n",
    "           bbox_to_anchor=(0.5, 0), loc='upper center', borderaxespad=3, frameon=False)\n",
    "##Save\n",
    "fileDir = './ExportFigures/'\n",
    "ipynbName = '220823_Multiomics-BMI-NatMed1stRevision_WHtR-LASSO-bcoef-ver2_'\n",
    "fileName = 'correlation-'+bwhtr+'-'+sex+'-vars_non-zero-in-all.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": [
    "### 5-6. Explained variance in WHtR by the variables retained in ≥1 model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> OLS linear regression model for significance: log(WHtR) ~ b0 + b1\\*Variable + b2\\*Sex + b3\\*BaseAge + b4\\*PCs  \n",
    "> Univariate model: log(WHtR) ~ b0 + b1\\*Variable  \n",
    ">\n",
    "> Of note, biological WHtR is included for comparison.  \n",
    "> As well as LASSO models, unstandardized WHtR is used for the dependent variable.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempD1 = {'Female':combiWHtR_F_bcoefs, 'Male':combiWHtR_M_bcoefs, 'BothSex':combiWHtR_B_bcoefs}\n",
    "tempD2 = {'Female':combiDF_F, 'Male':combiDF_M, 'BothSex':combiDF_B}\n",
    "tempD3 = {'Female':whtrDF_F, 'Male':whtrDF_M, 'BothSex':whtrDF_B}\n",
    "yvar = 'WHtR'\n",
    "bwhtr = 'CombiWHtR'\n",
    "topX = 30\n",
    "for sex in tempD1.keys():\n",
    "    print(sex+' model:')\n",
    "    \n",
    "    #Variables with non-zero beta-coefficient at least 1 of 10 models\n",
    "    tempDF = tempD1[sex]\n",
    "    tempL = tempDF.loc[tempDF['nZeros']==10].index.tolist()#To be excluded\n",
    "    \n",
    "    #Prepare DF for regressions\n",
    "    tempDF = tempD2[sex].drop(columns=tempL)\n",
    "    tempDF = tempDF.rename(columns={'log_Base'+bwhtr:bwhtr+' model'})\n",
    "    #Replace the standardized dependent variable with the unstandardized\n",
    "    tempDF1 = tempD3[sex]\n",
    "    tempDF['log_Base'+yvar] = tempDF1['log_Base'+yvar]\n",
    "    \n",
    "    #Independent variables for regressions (including bWHtR)\n",
    "    tempL = tempDF.drop(columns=whtrDF.columns.tolist()).columns.tolist()\n",
    "    \n",
    "    #Perform OLS linear regression\n",
    "    tempL1 = []#For sample size\n",
    "    tempL2 = []#For R2 in the univariate model\n",
    "    tempL3 = []#For R2\n",
    "    tempL4 = []#For beta-coefficient\n",
    "    tempL5 = []#For SE of beta-coefficient\n",
    "    tempL6 = []#For t-statistic\n",
    "    tempL7 = []#For residual degrees of freedom\n",
    "    tempL8 = []#For P-value\n",
    "    t_start = time.time()\n",
    "    for var in tempL:\n",
    "        #Rename independent variable\n",
    "        tempDF1 = tempDF.rename(columns={var: 'Variable'})\n",
    "        #Add a constant for the intercept -> Similar to R, smf automatically add a constant\n",
    "        \n",
    "        #Save sample size\n",
    "        tempL1.append(len(tempDF1))\n",
    "        \n",
    "        #Fit univariate model\n",
    "        formula = 'log_Base'+yvar+' ~ Variable'\n",
    "        fit_res = smf.ols(formula, data=tempDF1).fit()\n",
    "        #Save R2 [%]\n",
    "        tempL2.append(fit_res.rsquared*100)\n",
    "\n",
    "        #Fit full model\n",
    "        if sex=='BothSex':\n",
    "            formula = 'log_Base'+yvar+' ~ Variable + C(Sex) + BaseAge + PC1 + PC2 + PC3 + PC4 + PC5'\n",
    "        else:#Sex is invariant\n",
    "            formula = 'log_Base'+yvar+' ~ Variable + BaseAge + PC1 + PC2 + PC3 + PC4 + PC5'\n",
    "        fit_res = smf.ols(formula, data=tempDF1).fit()\n",
    "        #Save R2 [%]\n",
    "        tempL3.append(fit_res.rsquared*100)\n",
    "        #Save beta-coefficient of the variable\n",
    "        tempL4.append(fit_res.params['Variable'])\n",
    "        tempL5.append(fit_res.bse['Variable'])\n",
    "        #Save t-statistic of the variable\n",
    "        tempL6.append(fit_res.tvalues['Variable'])\n",
    "        #Save residual degrees of freedom\n",
    "        tempL7.append(int(fit_res.df_resid))\n",
    "        #Save P-value of the variable\n",
    "        tempL8.append(fit_res.pvalues['Variable'])\n",
    "    t_elapsed = time.time() - t_start\n",
    "    print('Elapsed time for', len(tempL), 'OLS linear regressions (including '+bwhtr+'):',\n",
    "          round(t_elapsed//60), 'min', round(t_elapsed%60, 1), 'sec')\n",
    "    \n",
    "    #Clean the results\n",
    "    tempDF = pd.DataFrame({'N':tempL1, 'UnivarR2':tempL2, 'R2':tempL3,\n",
    "                           'OLSbcoef':tempL4, 'OLSbcoefSE':tempL5,\n",
    "                           'tStat':tempL6, 'DoF':tempL7,  'Pval':tempL8},\n",
    "                          index=pd.Index(tempL, name='Variable'))\n",
    "    ##P-value adjustment by using Benjamini–Hochberg method\n",
    "    tempDF['AdjPval'] = multi.multipletests(tempDF['Pval'], alpha=0.05, method='fdr_bh',\n",
    "                                            is_sorted=False, returnsorted=False)[1]\n",
    "    ##Add the LASSO results\n",
    "    tempDF1 = tempD1[sex]\n",
    "    tempDF['LASSObcoef'] = tempDF1.loc[tempDF1.index.isin(tempL)]['Mean']#NaN in bWHtR reference\n",
    "    tempDF['LASSOnZeros'] = tempDF1.loc[tempDF1.index.isin(tempL)]['nZeros']#NaN in bWHtR reference\n",
    "    tempDF = tempDF.sort_values(by='AdjPval', ascending=True)\n",
    "    \n",
    "    #Save the cleaned DF\n",
    "    fileDir = './ExportData/'\n",
    "    ipynbName = '220823_Multiomics-BMI-NatMed1stRevision_WHtR-LASSO-bcoef-ver2_'\n",
    "    fileName = bwhtr+'-'+sex+'-OLS.tsv'\n",
    "    tempDF.to_csv(fileDir+ipynbName+fileName, sep='\\t', index=True)\n",
    "    \n",
    "    #Extact significant variables\n",
    "    tempDF = tempDF[tempDF['AdjPval']<0.05]\n",
    "    print('Variables significantly associated with WHtR (FDR < 0.05):', len(tempDF)-1)#bWHtR reference\n",
    "    ##Top X (+ bWHtR)\n",
    "    tempDF = tempDF.iloc[0:(topX+1)]\n",
    "    display(tempDF)\n",
    "    ##Category and color\n",
    "    tempL = []\n",
    "    for row_n in tempDF.index.tolist():\n",
    "        if row_n==bwhtr+' model':\n",
    "            tempL.append('Positive association')\n",
    "        else:\n",
    "            bcoef = tempDF.loc[row_n, 'OLSbcoef']\n",
    "            if bcoef>0:\n",
    "                tempL.append('Positive association')\n",
    "            elif bcoef<0:\n",
    "                tempL.append('Negative association')\n",
    "            else:\n",
    "                tempL.append('Error?')\n",
    "    tempDF['Association'] = tempL\n",
    "    tempD = {'Positive association':'tab:red', 'Negative association':'tab:blue'}\n",
    "    ##Plot\n",
    "    tempDF = tempDF.reset_index()\n",
    "    sns.set(style='ticks', font='Arial', context='talk')\n",
    "    plt.figure(figsize=(10, 4))\n",
    "    p = sns.barplot(data=tempDF, y='UnivarR2', x='Variable', hue='Association', hue_order=tempD.keys(),\n",
    "                    dodge=False, palette=tempD, edgecolor='k')\n",
    "    p.set(ylim=(0, 80), yticks=np.arange(0, 80, 15))#Fixed across omics\n",
    "    p.grid(axis='y', linestyle='--', color='black')\n",
    "    sns.despine()\n",
    "    plt.xticks(rotation=90, horizontalalignment='center')\n",
    "    plt.xlabel('')\n",
    "    plt.ylabel('Explained variance in WHtR [%]')\n",
    "    plt.legend(loc='upper right')\n",
    "    if sex=='BothSex':#Save\n",
    "        fileDir = './ExportFigures/'\n",
    "        ipynbName = '220823_Multiomics-BMI-NatMed1stRevision_WHtR-LASSO-bcoef-ver2_'\n",
    "        fileName = bwhtr+'-'+sex+'-OLS-top'+str(topX)+'.tif'\n",
    "        plt.gcf().savefig(fileDir+ipynbName+fileName, dpi=300, bbox_inches='tight', pad_inches=0.04,\n",
    "                          pil_kwargs={'compression':'tiff_lzw'})\n",
    "    plt.show()\n",
    "    print('')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6. Comparison b/w omics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6-1. Pearson's r"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 6-1-1. Actual values in LASSO models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Both sex cohort\n",
    "sex = 'BothSex'\n",
    "tempL1 = [metDF_B, protDF_B, chemDF_B, combiDF_B]\n",
    "tempL2 = ['log_BaseMetWHtR', 'log_BaseProtWHtR', 'log_BaseChemWHtR', 'log_BaseCombiWHtR']\n",
    "\n",
    "print('All variables')\n",
    "tempD = {'Metabolomics':'b', 'Proteomics':'r', 'Clinical labs':'g', 'Combined omics':'m'}\n",
    "tempDF = pd.DataFrame(columns=['Category', 'Variable1', 'Variable2', 'Pearson_r'])\n",
    "countL1 = []\n",
    "countL2 = []\n",
    "for df_i in range(len(tempL1)):\n",
    "    print(list(tempD.keys())[df_i])\n",
    "    #Compute correlation matrix\n",
    "    tempL = [item for sublist in [whtrDF.columns.tolist(), [tempL2[df_i]]] for item in sublist]\n",
    "    tempDF1 = tempL1[df_i].drop(columns=tempL)\n",
    "    tempDF1 = tempDF1.corr(method='pearson')\n",
    "    print(' - Combinations:', int(len(tempDF1)*(len(tempDF1)-1)/2))\n",
    "    #Extract lower triangle matrix\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",
    "    print(' - nrows:', len(tempDF1))\n",
    "    countL1.append(len(tempDF1))\n",
    "    print(' - |Pearson\\'s r| > 0.8:', len(tempDF1.loc[abs(tempDF1['Pearson_r'])>0.8]))\n",
    "    countL2.append(len(tempDF1.loc[abs(tempDF1['Pearson_r'])>0.8]))\n",
    "    #Clean\n",
    "    tempDF1['Category'] = list(tempD.keys())[df_i]\n",
    "    tempDF = pd.concat([tempDF, tempDF1], axis=0)\n",
    "print('Confirmation')\n",
    "display(tempDF['Category'].value_counts())\n",
    "\n",
    "#Distribution with violinplot\n",
    "sns.set(style='ticks', font='Arial', context='talk')\n",
    "plt.figure(figsize=(4, 3))\n",
    "ax = sns.violinplot(data=tempDF, x='Pearson_r', y='Category',\n",
    "                    order=list(tempD.keys()), palette=tempD, dodge=False, scale='width', inner='box')\n",
    "ax.set(xlim=(-1, 1), xticks=np.arange(-1, 1.1, 0.5))\n",
    "sns.despine()\n",
    "##Add annotation\n",
    "offset_x = 1.75\n",
    "ax.annotate('| Pearson\\'s '+r'$r$'+' | > 0.8', (offset_x, -1), fontsize='medium', annotation_clip=False,\n",
    "            verticalalignment='center', horizontalalignment='center')\n",
    "for df_i in range(len(tempD.keys())):\n",
    "    ax.annotate(f'{countL2[df_i]:,}'+' / '+f'{countL1[df_i]:,}'+' pairs',\n",
    "                (offset_x, df_i), fontsize='small', annotation_clip=False,\n",
    "                verticalalignment='center', horizontalalignment='center')\n",
    "plt.xlabel('Pearson\\'s '+r'$r$')\n",
    "plt.ylabel('')\n",
    "##Save\n",
    "fileDir = './ExportFigures/'\n",
    "ipynbName = '220823_Multiomics-BMI-NatMed1stRevision_WHtR-LASSO-bcoef-ver2_'\n",
    "fileName = 'correlation-'+sex+'-all-vars-pair.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": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Both sex cohort\n",
    "sex = 'BothSex'\n",
    "tempL1 = [metDF_B, protDF_B, chemDF_B, combiDF_B]\n",
    "tempL2 = [metWHtR_B_bcoefs, protWHtR_B_bcoefs, chemWHtR_B_bcoefs, combiWHtR_B_bcoefs]\n",
    "\n",
    "print('Variables with non-zero beta-coefficient in all 10 models')\n",
    "tempD = {'Metabolomics':'b', 'Proteomics':'r', 'Clinical labs':'g', 'Combined omics':'m'}\n",
    "tempDF = pd.DataFrame(columns=['Category', 'Variable1', 'Variable2', 'Pearson_r'])\n",
    "countL1 = []\n",
    "countL2 = []\n",
    "for df_i in range(len(tempL1)):\n",
    "    print(list(tempD.keys())[df_i])\n",
    "    ##Variables with non-zero beta-coefficient in all 10 models\n",
    "    tempDF1 = tempL2[df_i]\n",
    "    tempL = tempDF1.loc[tempDF1['nZeros']==0].index.tolist()\n",
    "    #Compute correlation matrix\n",
    "    tempDF1 = tempL1[df_i].loc[:, tempL]\n",
    "    tempDF1 = tempDF1.corr(method='pearson')\n",
    "    print(' - Combinations:', int(len(tempDF1)*(len(tempDF1)-1)/2))\n",
    "    #Extract lower triangle matrix\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",
    "    print(' - nrows:', len(tempDF1))\n",
    "    countL1.append(len(tempDF1))\n",
    "    print(' - |Pearson\\'s r| > 0.8:', len(tempDF1.loc[abs(tempDF1['Pearson_r'])>0.8]))\n",
    "    countL2.append(len(tempDF1.loc[abs(tempDF1['Pearson_r'])>0.8]))\n",
    "    #Clean\n",
    "    tempDF1['Category'] = list(tempD.keys())[df_i]\n",
    "    tempDF = pd.concat([tempDF, tempDF1], axis=0)\n",
    "print('Confirmation')\n",
    "display(tempDF['Category'].value_counts())\n",
    "\n",
    "#Distribution with violinplot\n",
    "sns.set(style='ticks', font='Arial', context='talk')\n",
    "plt.figure(figsize=(4, 3))\n",
    "ax = sns.violinplot(data=tempDF, x='Pearson_r', y='Category',\n",
    "                    order=list(tempD.keys()), palette=tempD, dodge=False, scale='width', inner='box')\n",
    "ax.set(xlim=(-1, 1), xticks=np.arange(-1, 1.1, 0.5))\n",
    "sns.despine()\n",
    "##Add annotation\n",
    "offset_x = 1.75\n",
    "ax.annotate('| Pearson\\'s '+r'$r$'+' | > 0.8', (offset_x, -1), fontsize='medium', annotation_clip=False,\n",
    "            verticalalignment='center', horizontalalignment='center')\n",
    "for df_i in range(len(tempD.keys())):\n",
    "    ax.annotate(f'{countL2[df_i]:,}'+' / '+f'{countL1[df_i]:,}'+' pairs',\n",
    "                (offset_x, df_i), fontsize='small', annotation_clip=False,\n",
    "                verticalalignment='center', horizontalalignment='center')\n",
    "plt.xlabel('Pearson\\'s '+r'$r$')\n",
    "plt.ylabel('')\n",
    "##Save\n",
    "fileDir = './ExportFigures/'\n",
    "ipynbName = '220823_Multiomics-BMI-NatMed1stRevision_WHtR-LASSO-bcoef-ver2_'\n",
    "fileName = 'correlation-'+sex+'-LASSO-vars-pair.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": [
    "### 6-2. Variables retained in all 10 models"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 6-2-1. vs. Metabolomics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempD1 = {'Female':metWHtR_F_bcoefs, 'Male':metWHtR_M_bcoefs, 'Both sex':metWHtR_B_bcoefs}\n",
    "tempD2 = {'Female':combiWHtR_F_bcoefs, 'Male':combiWHtR_M_bcoefs, 'Both sex':combiWHtR_B_bcoefs}\n",
    "tempT1 = ('MetWHtR', 'CombiWHtR')\n",
    "tempT2 = ('b', 'm')\n",
    "\n",
    "for sex in tempD1.keys():\n",
    "    print(sex+':')\n",
    "    tempDF1 = tempD1[sex]\n",
    "    tempDF2 = tempD2[sex]\n",
    "    \n",
    "    #Variables with non-zero beta-coefficient in all 10 models\n",
    "    tempS1 = set(tempDF1.loc[tempDF1['nZeros']==0].index)\n",
    "    tempS2 = set(tempDF2.loc[tempDF2['nZeros']==0].index)\n",
    "    \n",
    "    #Extract analytes in same omics type\n",
    "    tempS2 = tempS2 & set(tempDF1.index)\n",
    "    \n",
    "    #Common variables with non-zero beta-coefficient in all 10 models\n",
    "    tempS3 = tempS1 & tempS2\n",
    "    print('Common variables with non-zero beta-coefficient in all 10 models:\\n', tempS3)\n",
    "    \n",
    "    #Not common variables with non-zero beta-coefficient in all 10 models\n",
    "    tempS1 = tempS1 - tempS3\n",
    "    tempS2 = tempS2 - tempS3\n",
    "    \n",
    "    #The number of each subset\n",
    "    tempL = [len(tempS1), len(tempS2), len(tempS3)]\n",
    "    print('Each subset:', tempL)\n",
    "    \n",
    "    #Venn diagram\n",
    "    sns.set(font='Arial', context='talk')\n",
    "    venn2(subsets=tempL, set_labels=tempT1, set_colors=tempT2, alpha=0.7)\n",
    "    venn2_circles(subsets=tempL)\n",
    "    plt.title('—'+sex+' model—', fontdict={'fontsize':24})\n",
    "    plt.show()\n",
    "    print('')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 6-2-2. vs. Clinical labs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempD1 = {'Female':chemWHtR_F_bcoefs, 'Male':chemWHtR_M_bcoefs, 'Both sex':chemWHtR_B_bcoefs}\n",
    "tempD2 = {'Female':combiWHtR_F_bcoefs, 'Male':combiWHtR_M_bcoefs, 'Both sex':combiWHtR_B_bcoefs}\n",
    "tempT1 = ('ChemWHtR', 'CombiWHtR')\n",
    "tempT2 = ('g', 'm')\n",
    "\n",
    "for sex in tempD1.keys():\n",
    "    print(sex+':')\n",
    "    tempDF1 = tempD1[sex]\n",
    "    tempDF2 = tempD2[sex]\n",
    "    \n",
    "    #Variables with non-zero beta-coefficient in all 10 models\n",
    "    tempS1 = set(tempDF1.loc[tempDF1['nZeros']==0].index)\n",
    "    tempS2 = set(tempDF2.loc[tempDF2['nZeros']==0].index)\n",
    "    \n",
    "    #Extract analytes in same omics type\n",
    "    tempS2 = tempS2 & set(tempDF1.index)\n",
    "    \n",
    "    #Common variables with non-zero beta-coefficient in all 10 models\n",
    "    tempS3 = tempS1 & tempS2\n",
    "    print('Common variables with non-zero beta-coefficient in all 10 models:\\n', tempS3)\n",
    "    \n",
    "    #Not common variables with non-zero beta-coefficient in all 10 models\n",
    "    tempS1 = tempS1 - tempS3\n",
    "    tempS2 = tempS2 - tempS3\n",
    "    \n",
    "    #The number of each subset\n",
    "    tempL = [len(tempS1), len(tempS2), len(tempS3)]\n",
    "    print('Each subset:', tempL)\n",
    "    \n",
    "    #Venn diagram\n",
    "    sns.set(font='Arial', context='talk')\n",
    "    venn2(subsets=tempL, set_labels=tempT1, set_colors=tempT2, alpha=0.7)\n",
    "    venn2_circles(subsets=tempL)\n",
    "    plt.title('—'+sex+' model—', fontdict={'fontsize':24})\n",
    "    plt.show()\n",
    "    print('')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 6-2-3. vs. Proteomics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempD1 = {'Female':protWHtR_F_bcoefs, 'Male':protWHtR_M_bcoefs, 'Both sex':protWHtR_B_bcoefs}\n",
    "tempD2 = {'Female':combiWHtR_F_bcoefs, 'Male':combiWHtR_M_bcoefs, 'Both sex':combiWHtR_B_bcoefs}\n",
    "tempT1 = ('ProtWHtR', 'CombiWHtR')\n",
    "tempT2 = ('r', 'm')\n",
    "\n",
    "for sex in tempD1.keys():\n",
    "    print(sex+':')\n",
    "    tempDF1 = tempD1[sex]\n",
    "    tempDF2 = tempD2[sex]\n",
    "    \n",
    "    #Variables with non-zero beta-coefficient in all 10 models\n",
    "    tempS1 = set(tempDF1.loc[tempDF1['nZeros']==0].index)\n",
    "    tempS2 = set(tempDF2.loc[tempDF2['nZeros']==0].index)\n",
    "    \n",
    "    #Extract analytes in same omics type\n",
    "    tempS2 = tempS2 & set(tempDF1.index)\n",
    "    \n",
    "    #Common variables with non-zero beta-coefficient in all 10 models\n",
    "    tempS3 = tempS1 & tempS2\n",
    "    print('Common variables with non-zero beta-coefficient in all 10 models:\\n', tempS3)\n",
    "    \n",
    "    #Not common variables with non-zero beta-coefficient in all 10 models\n",
    "    tempS1 = tempS1 - tempS3\n",
    "    tempS2 = tempS2 - tempS3\n",
    "    \n",
    "    #The number of each subset\n",
    "    tempL = [len(tempS1), len(tempS2), len(tempS3)]\n",
    "    print('Each subset:', tempL)\n",
    "    \n",
    "    #Venn diagram\n",
    "    sns.set(font='Arial', context='talk')\n",
    "    venn2(subsets=tempL, set_labels=tempT1, set_colors=tempT2, alpha=0.7)\n",
    "    venn2_circles(subsets=tempL)\n",
    "    plt.title('—'+sex+' model—', fontdict={'fontsize':24})\n",
    "    plt.show()\n",
    "    print('')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# — End of this notebook —"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "arivale-py3 - Python",
   "language": "python",
   "name": "conda-env-arivale-py3-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  },
  "widgets": {
   "application/vnd.jupyter.widget-state+json": {
    "state": {},
    "version_major": 2,
    "version_minor": 0
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}