z_effect_on_shape.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "22d30a0f",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "HOME = os.environ[\"HOME\"]\n",
    "os.environ[\"CARDIAC_GWAS_REPO\"] = CARDIAC_GWAS_REPO = f\"{HOME}/01_repos/CardiacGWAS\"\n",
    "os.environ[\"CARDIAC_COMA_REPO\"] = CARDIAC_COMA_REPO = f\"{HOME}/01_repos/CardiacCOMA/\"\n",
    "os.environ[\"GWAS_REPO\"] = GWAS_REPO = f\"{HOME}/01_repos/GWAS_pipeline/\"\n",
    "\n",
    "MLRUNS_DIR = f\"{CARDIAC_COMA_REPO}/mlruns\"\n",
    "#os.chdir(CARDIAC_COMA_REPO)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "48550f22",
   "metadata": {},
   "outputs": [],
   "source": [
    "import mlflow\n",
    "from mlflow.tracking import MlflowClient\n",
    "\n",
    "import os, sys\n",
    "\n",
    "import torch\n",
    "import torch.nn.functional as F\n",
    "\n",
    "from CardiacCOMA.config.cli_args import overwrite_config_items\n",
    "from CardiacCOMA.config.load_config import load_yaml_config, to_dict\n",
    "from CardiacCOMA.utils.helpers import get_datamodule, get_lightning_module\n",
    "from CardiacCOMA.utils.mlflow_helpers import get_model_pretrained_weights\n",
    "from CardiacCOMA.utils.CardioMesh.CardiacMesh import transform_mesh\n",
    "\n",
    "import ipywidgets as widgets\n",
    "from ipywidgets import interact\n",
    "from IPython.display import Image, display, Markdown, clear_output\n",
    "\n",
    "import pandas as pd\n",
    "import shlex\n",
    "from subprocess import check_output\n",
    "\n",
    "import pickle as pkl\n",
    "import pytorch_lightning as pl\n",
    "\n",
    "from argparse import Namespace\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from IPython import embed\n",
    "sys.path.insert(0, '..')\n",
    "\n",
    "from copy import deepcopy\n",
    "from pprint import pprint\n",
    "\n",
    "from typing import List\n",
    "from tqdm import tqdm\n",
    "\n",
    "import pandas as pd\n",
    "\n",
    "import pyvista as pv\n",
    "from ipywidgets import interact, interactive, fixed, interact_manual\n",
    "\n",
    "from auxiliary import load_data\n",
    "from auxiliary import get_model_pretrained_weights\n",
    "import warnings\n",
    "\n",
    "from image_helpers import merge_pngs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "72955149",
   "metadata": {},
   "outputs": [],
   "source": [
    "def load_data(b):\n",
    "    global meshes, procrustes_transforms\n",
    "    print(\"Loading mesh data...\")\n",
    "    meshes = pkl.load(open(f\"{CARDIAC_COMA_REPO}/data/cardio/LV_meshes_at_ED_35k.pkl\", \"rb\"))\n",
    "    print(\"Mesh data loaded successfully.\")\n",
    "    \n",
    "    print(\"Loading Procrustes transforms...\")\n",
    "    procrustes_transforms = pkl.load(open(f\"{CARDIAC_COMA_REPO}/data/cardio/procrustes_transforms_35k.pkl\", \"rb\"))\n",
    "    print(\"Procrustes transform loaded successfully.\")\n",
    "    \n",
    "button = widgets.Button(\n",
    "    description='Load data',\n",
    "    disabled=False,\n",
    "    tooltip='Click me',\n",
    "    icon='check' # (FontAwesome names without the `fa-` prefix)\n",
    ")\n",
    "\n",
    "load_data(1)\n",
    "#button.on_click(load_data)\n",
    "#button    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fa07b9fe-75d6-4c5f-af9d-64c3984bc0f1",
   "metadata": {},
   "outputs": [],
   "source": [
    "ids = [ str(id) for id in meshes ] "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9322d562",
   "metadata": {},
   "outputs": [],
   "source": [
    "import random\n",
    "pv.set_plot_theme(\"document\")\n",
    "\n",
    "faces, _ = pkl.load(open(f\"{CARDIAC_COMA_REPO}/data/cardio/faces_and_downsampling_mtx_frac_0.1_LV.pkl\", \"rb\")).values()\n",
    "faces = np.c_[np.ones(faces.shape[0]) * 3, faces].astype(int)\n",
    "\n",
    "color_palette = list(pv.colors.color_names.values())\n",
    "random.shuffle(color_palette)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "163df1f3",
   "metadata": {},
   "source": [
    "# Select genetic locus"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "33ca1254",
   "metadata": {},
   "outputs": [],
   "source": [
    "good_runs_df = pd.read_csv(f\"{CARDIAC_GWAS_REPO}/results/good_runs.csv\")\n",
    "run_ids = good_runs_df.run_id.to_list()\n",
    "\n",
    "gwas_harmonized_pattern = \"data/other_gwas/preprocessed_files/{prefix}__{phenotype}.tsv\"\n",
    "gwas_selected_snps_pattern = \"data/other_gwas/preprocessed_files/{prefix}__{phenotype}__selected_snps.tsv\"\n",
    "COMA_GWAS_SUMMARY = \"results/gwas_loci_summary_across_runs.csv\"\n",
    "LOGP_PATH = \"results/log10p_for_selected_snps_across_gwas.csv\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6d4241d0",
   "metadata": {},
   "outputs": [],
   "source": [
    "good_runs_df = pd.read_csv(f\"{CARDIAC_GWAS_REPO}/results/good_runs.csv\")\n",
    "run_ids = good_runs_df.run_id.to_list()\n",
    "gwas_loci_summary_across_runs_df = pd.read_csv(COMA_GWAS_SUMMARY)\n",
    "\n",
    "# get index of best locus/variable\n",
    "idx = gwas_loci_summary_across_runs_df.groupby([\"region\"])[\"P\"].transform(min) == gwas_loci_summary_across_runs_df[\"P\"]\n",
    "\n",
    "best_association_per_region = gwas_loci_summary_across_runs_df[idx].sort_values(\"region\")\n",
    "best_snps = set(best_association_per_region.SNP)\n",
    "\n",
    "regions = { \n",
    "    f\"{assoc[1].region} ({assoc[1].P:.1e})\": assoc[1].region \n",
    "    for assoc in best_association_per_region.sort_values(\"P\").iterrows() \n",
    "} \n",
    "\n",
    "del regions[\"chr6_79 (4.5e-20)\"]\n",
    "\n",
    "_best_association_per_region = best_association_per_region.set_index(\"region\")\n",
    "# assoc = _best_association_per_region.loc[region]\n",
    "#run_id, z_variable = assoc.run, assoc.pheno[-4:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d9e804e5",
   "metadata": {},
   "outputs": [],
   "source": [
    "_best_association_per_region.head(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6419252c",
   "metadata": {},
   "outputs": [],
   "source": [
    "N_STEPS = 100\n",
    "STEP = 1 / N_STEPS\n",
    "quantiles = np.arange(N_STEPS+1)/N_STEPS\n",
    "\n",
    "quantile_ranges = [(quantiles[i], quantiles[i+1]) for i, q in enumerate(quantiles[:-1])]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a7526014-5c6b-4226-af86-f43a344b25c0",
   "metadata": {},
   "outputs": [],
   "source": [
    "#aligned_meshes = { id:  transform_mesh(mesh,  **procrustes_transforms[id]) for id, mesh in meshes.items() }\n",
    "#rmsd = {id: np.array([np.sqrt((p**2).sum()) for p in mesh]).mean() for id, mesh in aligned_meshes.items()}\n",
    "#pkl.dump(aligned_meshes, open(\"lved_aligned_meshes.pkl\", \"wb\"))\n",
    "#pkl.dump(rmsd, open(\"lved_rmsd.pkl\", \"wb\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ce91f794-e116-4927-b22d-8d1ad8ae160b",
   "metadata": {},
   "outputs": [],
   "source": [
    "aligned_meshes = pkl.load(open(\"lved_aligned_meshes.pkl\", \"rb\"))\n",
    "rmsd = pkl.load(open(\"lved_rmsd.pkl\", \"rb\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "437b4a4e-70e3-4e27-820d-9c51de553623",
   "metadata": {},
   "outputs": [],
   "source": [
    "# \n",
    "# scaled_meshes = []\n",
    "# msd = []\n",
    "#             \n",
    "# for id in ids:\n",
    "#     \n",
    "#     mesh = transform_mesh(\n",
    "#         meshes[id], \n",
    "#         **procrustes_transforms[id]\n",
    "#     )\n",
    "#         \n",
    "#     # mean-squared deviation\n",
    "#     msd.append(np.array([np.sqrt((p**2).sum()) for p in mesh]).mean())\n",
    "#     scaled_meshes.append(mesh)\n",
    "#     \n",
    "# scaled_meshes = np.array(scaled_meshes)\n",
    "# msd = np.array(msd)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d08243d5-5ba6-442a-8d4d-0357fe1e7130",
   "metadata": {},
   "outputs": [],
   "source": [
    "from typing import Sequence, Dict\n",
    "from z_effect_on_shape import avg_shape_in_q_range, get_ids_in_range, plot_mesh"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a0b16557-14d6-4510-bd25-3071e71ea314",
   "metadata": {},
   "outputs": [],
   "source": [
    "# TEST\n",
    "# ids_in_range = list(meshes.keys())[0:100]\n",
    "# avg_shape = np.array([meshes[id] / rmsd[id] for id in ids_in_range]).mean(0)\n",
    "# rmsd_in_range = np.array([rmsd[id] for id in ids_in_range]).mean()\n",
    "# avg_shape *= rmsd_in_range\n",
    "# avg_shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5c8ac039-3c86-4aa3-b8a2-4c57474e74c7",
   "metadata": {},
   "outputs": [],
   "source": [
    "exp_id = '1'\n",
    "\n",
    "regions = [\n",
    "   \"chr17_27\",\n",
    "   \"chr12_69\",\n",
    "   \"chr6_78\",\n",
    "   \"chr2_108\",\n",
    "   \"chr11_2\",\n",
    "   \"chr5_103\",\n",
    "   \"chr1_124\",\n",
    "   \"chr6_84\",\n",
    "   \"chr2_69\",\n",
    "   \"chr12_67\"\n",
    "]\n",
    "\n",
    "q_ranges = [(0.00, 0.01), (0.095, 0.105), (0.45, 0.55), (0.895, 0.905), (0.99, 1.0)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "68bb7731-0d89-44b0-9ae7-5ca6abdb005b",
   "metadata": {},
   "outputs": [],
   "source": [
    "VERBOSE = False\n",
    "WARNINGS = False \n",
    "\n",
    "if WARNINGS:\n",
    "    warnings.filterwarnings('default')    \n",
    "else:\n",
    "    warnings.filterwarnings('ignore')\n",
    "\n",
    "    \n",
    "for run_id in tqdm(sorted(os.listdir(f\"{CARDIAC_COMA_REPO}/mlruns/1\"))):\n",
    "    \n",
    "    try:\n",
    "        z_filepath = f\"{MLRUNS_DIR}/{exp_id}/{run_id}/artifacts/output/latent_vector.csv\"\n",
    "        zs = pd.read_csv(z_filepath).set_index(\"ID\").columns\n",
    "    except:\n",
    "        continue\n",
    "        \n",
    "    for z in zs:\n",
    "                \n",
    "        #assoc = _best_association_per_region.loc[region]\n",
    "        #run_id, z = assoc.run, assoc.pheno[-4:]\n",
    "    \n",
    "        z_effect_figs_dir = f\"{CARDIAC_COMA_REPO}/mlruns/1/{run_id}/artifacts/z_effect_on_shape\"\n",
    "        os.makedirs(z_effect_figs_dir, exist_ok=True)\n",
    "        \n",
    "        filenames = { (q0, q1): f\"{z_effect_figs_dir}/{z}_{q0}-{q1}.png\" for q0, q1 in q_ranges }\n",
    "    \n",
    "        for (q0, q1), filename in filenames.items():\n",
    "              \n",
    "            if os.path.exists(filename):\n",
    "                if VERBOSE: print(f\"File {filename} already exists.\")\n",
    "                continue                        \n",
    "            \n",
    "            avg_mesh = avg_shape_in_q_range(meshes=aligned_meshes, rmsd=rmsd, quantile_range=(q0, q1), exp_id=exp_id, run_id=run_id, z=z)  \n",
    "            \n",
    "            if np.isnan(avg_mesh).all(): break\n",
    "                                \n",
    "            ofilename = filenames[(q0, q1)]\n",
    "            plot_mesh(avg_mesh, faces, ofilename)                                    \n",
    "           \n",
    "        try:            \n",
    "            filename = f\"{z_effect_figs_dir}/{z}.png\"\n",
    "            if not os.path.exists(filename):\n",
    "                merge_pngs(filenames.values(), output_png=filename, how=\"horizontally\")\n",
    "        except:\n",
    "            pass\n",
    "        \n",
    "    if VERBOSE:\n",
    "        print(f\"{run_id[:5]}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "70a3d3c0-4a9c-4fcb-9176-2778181ecae6",
   "metadata": {},
   "outputs": [],
   "source": [
    "for region in tqdm(regions):\n",
    "  \n",
    "  assoc = _best_association_per_region.loc[region]\n",
    "  run_id, z = assoc.run, assoc.pheno[-4:]\n",
    "\n",
    "  df = pd.read_csv(f\"{MLRUNS_DIR}/{exp_id}/{run_id}/artifacts/output/latent_vector.csv\").set_index(\"ID\")\n",
    "\n",
    "  filenames = []\n",
    "\n",
    "  for q0, q1 in q_ranges:  \n",
    "    \n",
    "    filename = f\"results/figs/{region}_{q0}-{q1}.png\"    \n",
    "    filenames.append(filename)  \n",
    "    \n",
    "    if os.path.exists(filename):\n",
    "        print(f\"File {filename} already exists.\")\n",
    "        continue\n",
    "    \n",
    "    z_bounds = df.quantile([q0, q1])[z]                \n",
    "        \n",
    "    ids = list(df[ \n",
    "        (z_bounds[q0] < df[z]) & (df[z] < z_bounds[q1])\n",
    "    ][z].index) \n",
    "      \n",
    "    scaled_meshes = []\n",
    "    msd = []\n",
    "    \n",
    "    for id in ids:\n",
    "        mesh = transform_mesh(\n",
    "            meshes[str(id)], \n",
    "            **procrustes_transforms[str(id)]\n",
    "        )\n",
    "        \n",
    "        # mean-squared deviation\n",
    "        msd.append(np.array([np.sqrt((p**2).sum()) for p in mesh]).mean())\n",
    "        scaled_meshes.append(mesh)\n",
    "    \n",
    "    scaled_meshes = np.array(scaled_meshes)\n",
    "    msd = np.array(msd)\n",
    "            \n",
    "    avg_mesh = np.array([scaled_meshes[i] / msd[i] for i, _ in enumerate(scaled_meshes)]).mean(0)\n",
    "    avg_mesh = avg_mesh * msd.mean()\n",
    "                \n",
    "    pv.set_plot_theme(\"document\")\n",
    "    pl = pv.Plotter(off_screen=True, notebook=False)\n",
    "    \n",
    "    pl.camera.position = (300, 0.0, 0.0)\n",
    "    pl.camera.azimuth = 95\n",
    "    \n",
    "    mesh = pv.PolyData(avg_mesh, faces)\n",
    "    pl.add_mesh(mesh, show_edges=False, point_size=1.5, color=color_palette[0], opacity=0.5)\n",
    "    \n",
    "    pl.screenshot(filename);\n",
    "    \n",
    "  \n",
    "  merge_pngs(filenames, output_png=f\"{region}.png\", how=\"horizontally\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8d4edff2",
   "metadata": {},
   "outputs": [],
   "source": [
    "exp_id = '1'\n",
    "# region = \"chr17_27\"\n",
    "region = \"chr12_69\"\n",
    "assoc = _best_association_per_region.loc[region]\n",
    "run_id, z = assoc.run, assoc.pheno[-4:]\n",
    "filename = \"latent_vector.csv\"\n",
    "df = pd.read_csv(f\"{MLRUNS_DIR}/{exp_id}/{run_id}/artifacts/output/{filename}\")\n",
    "df = df.set_index(\"ID\")\n",
    "\n",
    "for q0, q1 in quantile_ranges:\n",
    "  \n",
    "    z_bounds = df.quantile([q0, q1])[z]                \n",
    "    \n",
    "    ids = list(df[ \n",
    "        (z_bounds[q0] < df[z]) & (df[z] < z_bounds[q1])\n",
    "    ][z].index) \n",
    "    \n",
    "    scaled_meshes = []\n",
    "    msd = []\n",
    "    \n",
    "    for id in ids:\n",
    "        mesh = transform_mesh(\n",
    "            meshes[str(id)], \n",
    "            **procrustes_transforms[str(id)]\n",
    "        )\n",
    "        \n",
    "        # mean-squared deviation\n",
    "        msd.append(np.array([np.sqrt((p**2).sum()) for p in mesh]).mean())\n",
    "        scaled_meshes.append(mesh)\n",
    "    \n",
    "    scaled_meshes = np.array(scaled_meshes)\n",
    "    \n",
    "    msd = np.array(msd)            \n",
    "    avg_mesh = np.array([scaled_meshes[i] / msd[i] for i, _ in enumerate(scaled_meshes)]).mean(0)\n",
    "    avg_mesh = avg_mesh * msd.mean()\n",
    "    \n",
    "    pl = pv.Plotter(notebook=True, off_screen=False, polygon_smoothing=False)\n",
    "    mesh = pv.PolyData(avg_mesh, faces)\n",
    "    pl.add_mesh(mesh, show_edges=False, point_size=1.5, color=color_palette[0], opacity=0.5)\n",
    "    \n",
    "    clear_output()\n",
    "    print()\n",
    "    pl.show(interactive=True, interactive_update=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fb17bc37",
   "metadata": {},
   "outputs": [],
   "source": [
    "exp_id = '1'\n",
    "STEP = 0.01\n",
    "\n",
    "region_w = widgets.Select(options=regions, description=\"Locus: \\n\",)\n",
    "display(region_w)\n",
    "region = region_w.value\n",
    "\n",
    "# _best_association_per_region = best_association_per_region.set_index(\"region\")\n",
    "assoc = _best_association_per_region.loc[region]\n",
    "run_id, z = assoc.run, assoc.pheno[-4:]\n",
    "\n",
    "quantile_inf_w=widgets.FloatSlider(min=0.00, max=0.99, step=STEP)\n",
    "display(quantile_inf_w)    \n",
    "\n",
    "button = widgets.Button(description=\"Plot mesh\")\n",
    "out = widgets.Output()\n",
    "\n",
    "def plot_mesh(b):\n",
    "    \n",
    "    filename = \"latent_vector.csv\"\n",
    "    df = pd.read_csv(f\"{MLRUNS_DIR}/{exp_id}/{run_id}/artifacts/output/{filename}\")\n",
    "    df = df.set_index(\"ID\")\n",
    "    \n",
    "    # z_w = widgets.Select(options=[f\"z{str(i).zfill(3)}\" for i in range(df.shape[1])])    \n",
    "    \n",
    "    with out:\n",
    "\n",
    "        quantile_range = quantile_inf_w.value\n",
    "        z_bounds = df.quantile([quantile_range, STEP+quantile_range])[z]                \n",
    "        \n",
    "        ids = list(df[ \n",
    "            (z_bounds[quantile_range] < df[z]) & (df[z] < z_bounds[quantile_range+STEP])\n",
    "        ][z].index)    \n",
    "                \n",
    "        # COMPUTE AVERAGE MESH                        \n",
    "        avg_mesh = np.array([\n",
    "            transform_mesh(\n",
    "                meshes[str(id)], \n",
    "                **procrustes_transforms[str(id)]\n",
    "            ) for id in ids\n",
    "        ]).mean(axis=0)\n",
    "\n",
    "        clear_output()        \n",
    "        \n",
    "        pl = pv.Plotter(notebook=True, off_screen=False, polygon_smoothing=False)\n",
    "        mesh = pv.PolyData(avg_mesh, faces)\n",
    "        pl.add_mesh(mesh, show_edges=False, point_size=1.5, color=color_palette[0], opacity=0.5)\n",
    "        pl.show(interactive=True, interactive_update=True)\n",
    "        \n",
    "button.on_click(plot_mesh)\n",
    "# displaying button and its output together\n",
    "widgets.VBox([button,out])    "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "184eca27",
   "metadata": {},
   "source": [
    "___"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f87f85eb",
   "metadata": {},
   "source": [
    "# Select `run_id` / `z` variable"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1be5bd8f",
   "metadata": {},
   "outputs": [],
   "source": [
    "from PIL import Image\n",
    "\n",
    "exp_id = '1'\n",
    "STEP=0.01\n",
    "run_id_w = widgets.Select(options=sorted(run_ids))\n",
    "display(run_id_w)\n",
    "\n",
    "filename = \"latent_vector.csv\"\n",
    "df = pd.read_csv(f\"{MLRUNS_DIR}/{exp_id}/{run_id}/artifacts/output/{filename}\")\n",
    "df = df.set_index(\"ID\")\n",
    "\n",
    "z_w = widgets.SelectionSlider(options=[f\"z{str(i).zfill(3)}\" for i in range(df.shape[1])])\n",
    "display(z_w)\n",
    "\n",
    "quantile_range_w=widgets.FloatSlider(min=0, max=0.99, step=STEP)\n",
    "display(quantile_range_w)\n",
    "\n",
    "button = widgets.Button(description=\"Plot mesh\")\n",
    "\n",
    "out = widgets.Output()\n",
    "\n",
    "\n",
    "def plot_mesh(b):\n",
    "    \n",
    "    with out:\n",
    "\n",
    "        z = z_w.value\n",
    "        run_id = run_id_w.value\n",
    "        quantile_range = quantile_range_w.value\n",
    "        z_bounds = df.quantile([quantile_range_w.value, quantile_range_w.value+STEP])[z]                \n",
    "        \n",
    "        ids = list(df[ \n",
    "            (z_bounds[quantile_range_w.value] < df[z]) & (df[z] < z_bounds[quantile_range_w.value+STEP])\n",
    "        ][z].index) \n",
    "        \n",
    "        print(len(ids))\n",
    "                \n",
    "        clear_output()      \n",
    "        print(z, run_id)\n",
    "        manhattan_file = f\"{CARDIAC_COMA_REPO}/mlruns/1/{run_id}/artifacts/GWAS_adj_10PCs/figures/GWAS__{z}__1_{run_id}__manhattan.png\"\n",
    "        qq_file = f\"{CARDIAC_COMA_REPO}/mlruns/1/{run_id}/artifacts/GWAS_adj_10PCs/figures/GWAS__{z}__1_{run_id}__QQ-plot.png\"\n",
    "        \n",
    "        # COMPUTE AVERAGE MESH                        \n",
    "        avg_mesh = np.array([\n",
    "            transform_mesh(\n",
    "                meshes[str(id)], \n",
    "                **procrustes_transforms[str(id)]\n",
    "            ) for id in ids\n",
    "        ]).mean(axis=0)\n",
    "\n",
    "                      \n",
    "        pl = pv.Plotter(notebook=True, off_screen=False, polygon_smoothing=False)\n",
    "        mesh = pv.PolyData(avg_mesh, faces)\n",
    "        pl.add_mesh(mesh, show_edges=False, point_size=1.5, color=color_palette[0], opacity=0.5)\n",
    "        pv.camera\n",
    "        pl.show(interactive=True, interactive_update=True)\n",
    "        \n",
    "        display(Image.open(manhattan_file))\n",
    "        # display(Image.open(qq_file))\n",
    "        \n",
    "button.on_click(plot_mesh)\n",
    "# displaying button and its output together\n",
    "widgets.VBox([button,out])        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3f0bfd4d",
   "metadata": {},
   "outputs": [],
   "source": [
    "pl = pv.Plotter(notebook=True, off_screen=False, polygon_smoothing=False)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}