main_inference_heatmaps.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This tutorial focuses on how to make MSI predictions at both slide-level and tile-level with the trained MSI-SEER models. We provide the model weights trained using Yonsei-1 for colorectal cancer slides and using the combined data (TCGA-STAD and Yonsei-Classic) for gastric cancer. The Python script file \"script_msi_seer_train_inference.py\" contains the codes that train the MSI-SEER model using the image features computed by the pre-trained deep learning models. We have included two example slides that were used to generate Figure 4 (the tile-level prediction heatmaps) in this repository so that readers can understand how our model works on whole slide images and check the reproducibility of our model. Please note that MSI-SEER is based on Monte Carlo sampling to implement both training and testing procedures, so there may be some randomness in prediction results. However, prediction results should not be too different from those reported in the paper."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [],
   "source": [
    "#- load the packages \n",
    "import numpy as np\n",
    "import os\n",
    "\n",
    "from models.wsl_binary_classifier_fea_vi_agg_ensemble import wsl_classifier as dgp_rf_agg_ens\n",
    "from models.wsl_binary_classifier_fea_vi_agg_ensemble import cal_unc_quntities\n",
    "\n",
    "from models.utils import get_MSI_data_with_tileinfo, data_feature_mat, seed_everything\n",
    "\n",
    "import copy\n",
    "\n",
    "#- fix the random seed\n",
    "seed_everything(11111)\n",
    "\n",
    "N_MSIDETECT_models = 9    \n",
    "\n",
    "#- options for MSI-SEER\n",
    "from types import SimpleNamespace   \n",
    "setting = {'batch_size': 8,\n",
    "    'max_epoch': 100,\n",
    "    'minimum_epoch': 75,\n",
    "    'sub_Ni': 300,\n",
    "    'sel_instances': 'sampling',\n",
    "    'nMCsamples': 10 ,   \n",
    "    'alpha': 0.5,\n",
    "    'n_layers': 6,\n",
    "    'n_RFs': 100,\n",
    "    'ker_type': 'arccosin',\n",
    "    'exp_ensemble': True,\n",
    "    'n_runs': 10,\n",
    "    'flag_mean_function': True, \n",
    "    'iter_print': True,\n",
    "    'save_directory': './model_weights',\n",
    "    'w_pos': 1.0,\n",
    "    'model_ref_path': './model_weights/prior_means_fromMSIDETECT/'\n",
    "}\n",
    "setting = SimpleNamespace(**setting)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#- Codes for making both tile-level and tile-level predictions on the selected examples\n",
    "PATH_FEATURES = './sample_wsi_features'\n",
    "\n",
    "str_save_path = './results/predictions/examples'\n",
    "os.makedirs(str_save_path, exist_ok=True)\n",
    "\n",
    "str_trn_data = 'GC_Combined'\n",
    "\n",
    "Yest_list = [] \n",
    "for ith_run in range(N_MSIDETECT_models):\n",
    "    str_model = 'MODEL_' + str(ith_run)  \n",
    "\n",
    "    str_save_each_model_path =  os.path.join(str_save_path, str_model)\n",
    "    if not os.path.exists(str_save_each_model_path):\n",
    "        os.makedirs(str_save_each_model_path)   \n",
    "\n",
    "    # data (image features) loading\n",
    "    img_names_tst, img_tilenames_tst, X_tst, _, _ \\\n",
    "            = get_MSI_data_with_tileinfo(os.path.join(PATH_FEATURES, str_model))                                    \n",
    "        \n",
    "    if ith_run == 0:\n",
    "        img_names_tst_ref = copy.deepcopy(img_names_tst)\n",
    "    else:   \n",
    "        idx_matched = [img_names_tst.index(elm) for elm in img_names_tst_ref]\n",
    "        X_tst = [X_tst[idx] for idx in idx_matched]\n",
    "        img_tilenames_tst = [img_tilenames_tst[idx] for idx in idx_matched]\n",
    "    \n",
    "    N_total = len(img_names_tst_ref)\n",
    "         \n",
    "    X_tst_cs = data_feature_mat \\\n",
    "    (sample_ids=img_names_tst_ref, tile_names=img_tilenames_tst, data_mat=X_tst)\n",
    "\n",
    "    # create the MSI model   \n",
    "    str_trn_model_name = str_trn_data + '/MODEL_' + str(ith_run)\n",
    "          \n",
    "    model_sepbest = dgp_rf_agg_ens(X_tst_cs, None, setting, str_trndata=str_trn_model_name)\n",
    "    \n",
    "    #- predictions \n",
    "    Ytst_probs_V1, _ = model_sepbest.predict\\\n",
    "                    (np.arange(len(img_names_tst)), data_set_=X_tst_cs, save_pred_path=str_save_each_model_path)\n",
    "                    \n",
    "    Yest_list.append(Ytst_probs_V1)\n",
    "    \n",
    "# inference in ensemble learning\n",
    "Yest_all = []\n",
    "for mn_sub in range(N_total):\n",
    "    Yest_cur = []\n",
    "    for ith_run in range(N_MSIDETECT_models):  # N_MSIDETECT_models\n",
    "        Yest_cur.append(Yest_list[ith_run][mn_sub])\n",
    "\n",
    "    Yest_all.append(np.hstack(Yest_cur))\n",
    "\n",
    "# calculate the aggregated MSI prediction probabilities and Bayesian confidence scores \n",
    "[Ytst_mean, unc_aleat, unc_epist] = cal_unc_quntities(Yest_all)\n",
    "BCS = 1 - 2*np.sqrt(unc_aleat + unc_epist)               "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import os\n",
    "from PIL import Image \n",
    "\n",
    "N_MSIDETECT_models = 9  \n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import cblind as cb\n",
    "\n",
    "from models.utils import fill_heatmap_grid_upscale\n",
    "from scipy.ndimage.filters import gaussian_filter\n",
    "\n",
    "img_list = ['Item_1', 'Item_2']\n",
    "\n",
    "str_save_path = './results/predictions/examples/'\n",
    "\n",
    "sigmoid_np = lambda x: 1.0 / (1.0 + np.exp(-x))\n",
    "\n",
    "for img in img_list: \n",
    "    #- calculate the MSI-prediction probabilities at tile-level               \n",
    "    for ith_run in range(N_MSIDETECT_models):\n",
    "        resolved_mat_path = os.path.join(str_save_path, f'MODEL_{ith_run}',  img + '_patch_maps.npz')\n",
    "\n",
    "        patch_mat = np.load(resolved_mat_path, allow_pickle=True)\n",
    "\n",
    "        if ith_run ==0:\n",
    "            patch_ids = np.reshape(patch_mat['arr_0'], [-1, 1])\n",
    "            patch_Fout = patch_mat['arr_1']\n",
    "        else:\n",
    "            idx_matched = np.reshape(np.array([np.where(patch_mat['arr_0']==elm)[0] for elm in patch_ids]), [-1])\n",
    "            assert np.sum(idx_matched == np.array(range(len(idx_matched)))) == len(idx_matched)\n",
    "\n",
    "            patch_Fout = np.concatenate((patch_Fout,  patch_mat['arr_1']), axis=1)\n",
    "    \n",
    "    patch_preds = np.mean(sigmoid_np(patch_Fout), axis=1, keepdims=True) \n",
    "    \n",
    "    # Get WSI details\n",
    "    wsi_width, wsi_height = np.load(f'./sample_wsi_features/{img}_wsiinfo.npy')\n",
    "\n",
    "    # Get thumbnail\n",
    "    thumbnail = Image.open(f'./sample_wsi_features/{img}_thumb.png')\n",
    "   \n",
    "    #- generate the heatmap  \n",
    "    n = int(patch_ids[0][0].split(\"_dy\")[-1].split(patch_ids[0][0][-4:])[0])            \n",
    "                    \n",
    "    width, height = thumbnail.size[0], thumbnail.size[1]             \n",
    "\n",
    "    grid_dims = (np.array((wsi_width, wsi_height)) / n).astype(np.int32)\n",
    "\n",
    "    heatmap_image = fill_heatmap_grid_upscale(\\\n",
    "        [wsi_width, wsi_height], n, [width, height], patch_ids, patch_preds)                                           \n",
    "\n",
    "    heatmap_image_np = np.array(heatmap_image)                                             \n",
    "\n",
    "    #- tile-level prediction heatmap\n",
    "    plt.figure()\n",
    "\n",
    "    plt.imshow(thumbnail)\n",
    "            \n",
    "    scale_factor_w = np.sqrt(width/grid_dims[1])\n",
    "    filtered_arr = gaussian_filter(heatmap_image_np, sigma=scale_factor_w)\n",
    "            \n",
    "    plt.imshow(np.ma.array(filtered_arr, mask=heatmap_image_np==0), \\\n",
    "        interpolation='hanning', alpha=0.4, cmap=cb.cbmap(\"cb.solstice\", nbin=20))\n",
    "            \n",
    "    cbar = plt.colorbar()\n",
    "    cbar.set_label(\"MSI-H predictive probability\")\n",
    "\n",
    "    plt.tight_layout(w_pad=0)\n",
    "    plt.axis('off')\n",
    "\n",
    "    plt.show()\n",
    "\n",
    "    #- histogram of the MSI prediction probabilities in the slide\n",
    "    fig, ax = plt.subplots()\n",
    "\n",
    "    ax.hist(patch_preds, density=False, bins=50, alpha=0.75)\n",
    "\n",
    "    plt.xlabel(\"MSI-H predictive probability\")\n",
    "    plt.show()\n",
    "    "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "msi_seer",
   "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.19"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}