7.figures/SupplementaryFigure_FeatureImportance.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "07795c5d-b696-4984-b5d4-4573df5189c2",
   "metadata": {},
   "source": [
    "## Visualizing feature importance in multiclass models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "a76f063a-d371-42b6-bfea-ae57f5bae19b",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning message:\n",
      "“package ‘ggplot2’ was built under R version 4.2.3”\n",
      "Warning message:\n",
      "“package ‘patchwork’ was built under R version 4.2.3”\n"
     ]
    }
   ],
   "source": [
    "suppressPackageStartupMessages(library(dplyr))\n",
    "suppressPackageStartupMessages(library(ggplot2))\n",
    "suppressPackageStartupMessages(library(patchwork))\n",
    "suppressPackageStartupMessages(library(ComplexHeatmap))\n",
    "\n",
    "# Load variables important for plotting (e.g., themes, phenotypes, etc.)\n",
    "source(\"themes.r\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "cde2ac8d-e357-4e9a-8a68-f94df0c2dd83",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set output directory\n",
    "output_dir <- file.path(\"figures\", \"individual_coefficient_heatmaps\")\n",
    "\n",
    "heatmap_real_file <- file.path(output_dir, \"compiled_real_coefficient_heatmaps.png\")\n",
    "heatmap_shuffled_file <- file.path(output_dir, \"compiled_shuffled_coefficient_heatmaps.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "30db84d4-e37a-4746-b4b1-90c5df36aeb2",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\u001b[1m\u001b[22mNew names:\n",
      "\u001b[36m•\u001b[39m `` -> `...1`\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] 86220     5\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<table class=\"dataframe\">\n",
       "<caption>A tibble: 3 × 5</caption>\n",
       "<thead>\n",
       "\t<tr><th scope=col>Feature_Name</th><th scope=col>Phenotypic_Class</th><th scope=col>Coefficent_Value</th><th scope=col>shuffled</th><th scope=col>feature_type</th></tr>\n",
       "\t<tr><th scope=col>&lt;chr&gt;</th><th scope=col>&lt;chr&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;chr&gt;</th><th scope=col>&lt;chr&gt;</th></tr>\n",
       "</thead>\n",
       "<tbody>\n",
       "\t<tr><td>CP__AreaShape_Area</td><td>ADCCM    </td><td> 0.00822653</td><td>False</td><td>CP</td></tr>\n",
       "\t<tr><td>CP__AreaShape_Area</td><td>Anaphase </td><td>-0.01862480</td><td>False</td><td>CP</td></tr>\n",
       "\t<tr><td>CP__AreaShape_Area</td><td>Apoptosis</td><td>-0.02102483</td><td>False</td><td>CP</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "A tibble: 3 × 5\n",
       "\\begin{tabular}{lllll}\n",
       " Feature\\_Name & Phenotypic\\_Class & Coefficent\\_Value & shuffled & feature\\_type\\\\\n",
       " <chr> & <chr> & <dbl> & <chr> & <chr>\\\\\n",
       "\\hline\n",
       "\t CP\\_\\_AreaShape\\_Area & ADCCM     &  0.00822653 & False & CP\\\\\n",
       "\t CP\\_\\_AreaShape\\_Area & Anaphase  & -0.01862480 & False & CP\\\\\n",
       "\t CP\\_\\_AreaShape\\_Area & Apoptosis & -0.02102483 & False & CP\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "A tibble: 3 × 5\n",
       "\n",
       "| Feature_Name &lt;chr&gt; | Phenotypic_Class &lt;chr&gt; | Coefficent_Value &lt;dbl&gt; | shuffled &lt;chr&gt; | feature_type &lt;chr&gt; |\n",
       "|---|---|---|---|---|\n",
       "| CP__AreaShape_Area | ADCCM     |  0.00822653 | False | CP |\n",
       "| CP__AreaShape_Area | Anaphase  | -0.01862480 | False | CP |\n",
       "| CP__AreaShape_Area | Apoptosis | -0.02102483 | False | CP |\n",
       "\n"
      ],
      "text/plain": [
       "  Feature_Name       Phenotypic_Class Coefficent_Value shuffled feature_type\n",
       "1 CP__AreaShape_Area ADCCM             0.00822653      False    CP          \n",
       "2 CP__AreaShape_Area Anaphase         -0.01862480      False    CP          \n",
       "3 CP__AreaShape_Area Apoptosis        -0.02102483      False    CP          "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Load model coefficients\n",
    "coef_dir <- file.path(\"../4.interpret_model/coefficients/\")\n",
    "coef_file <- file.path(coef_dir, \"compiled_coefficients.tsv\")\n",
    "\n",
    "coef_df <- readr::read_tsv(\n",
    "    coef_file,\n",
    "    col_types = readr::cols(\n",
    "        .default = \"d\",\n",
    "        \"Feature_Name\" = \"c\",\n",
    "        \"Phenotypic_Class\" = \"c\",\n",
    "        \"feature_type\" = \"c\",\n",
    "        \"shuffled\" = \"c\"\n",
    "    )\n",
    ") %>%\n",
    "    dplyr::select(!`...1`) \n",
    "\n",
    "print(dim(coef_df))\n",
    "head(coef_df, 3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "4074aced-f249-4029-bf51-8da9e5103a44",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set constants for how to split feature name into interpretable parts\n",
    "first_component_split <- c(\n",
    "    \"feature_space\",\n",
    "    \"feature_name\"\n",
    ")\n",
    "\n",
    "dp_feature_components <- c(\n",
    "    \"model\",\n",
    "    \"component_idx\"\n",
    ")\n",
    "\n",
    "cp_feature_components <- c(\n",
    "    \"feature_group\",\n",
    "    \"measurement\",\n",
    "    \"channel\", \n",
    "    \"parameter1\", \n",
    "    \"parameter2\",\n",
    "    \"parameter3\"\n",
    ")\n",
    "\n",
    "metadata_annotation_options <- list(\n",
    "    \"CP\" = cp_feature_components,\n",
    "    \"DP\" = dp_feature_components,\n",
    "    \"CP_and_DP\" = dp_feature_components\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "dd973800-e19a-489d-9b0b-029ff143d74e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] \"Transforming model coefficients: CP__False\"\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning message:\n",
      "“\u001b[1m\u001b[22mExpected 6 pieces. Missing pieces filled with `NA` in 105 rows [1, 2, 3, 4, 5,\n",
      "6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, ...].”\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] \"Transforming model coefficients: CP_and_DP__False\"\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning message:\n",
      "“\u001b[1m\u001b[22mExpected 2 pieces. Additional pieces discarded in 138 rows [3, 4, 5, 6, 7, 8,\n",
      "26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, ...].”\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] \"Transforming model coefficients: DP__False\"\n",
      "[1] \"Transforming model coefficients: CP__True\"\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning message:\n",
      "“\u001b[1m\u001b[22mExpected 6 pieces. Missing pieces filled with `NA` in 105 rows [1, 2, 3, 4, 5,\n",
      "6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, ...].”\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] \"Transforming model coefficients: CP_and_DP__True\"\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning message:\n",
      "“\u001b[1m\u001b[22mExpected 2 pieces. Additional pieces discarded in 138 rows [3, 4, 5, 6, 7, 8,\n",
      "26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, ...].”\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] \"Transforming model coefficients: DP__True\"\n"
     ]
    }
   ],
   "source": [
    "# Create matrices for each combination of shuffled and feature type\n",
    "feature_matrix_list <- list()\n",
    "for (shuffled_value in unique(coef_df$shuffled)) {\n",
    "    for (feature_type_value in unique(coef_df$feature_type)) {\n",
    "        list_name_index <- paste0(feature_type_value, \"__\", shuffled_value)\n",
    "        print(paste(\"Transforming model coefficients:\", list_name_index))\n",
    "        \n",
    "        # Subset the coeficient dataframe\n",
    "        coef_subset <- coef_df %>%\n",
    "            dplyr::filter(\n",
    "                shuffled == !!shuffled_value,\n",
    "                feature_type == !!feature_type_value\n",
    "            ) %>%\n",
    "            dplyr::select(-c(\"shuffled\", \"feature_type\")) %>%\n",
    "            tidyr::pivot_wider(names_from = Feature_Name, values_from = Coefficent_Value)\n",
    "\n",
    "        # Make data ready for matrix transformation\n",
    "        coef_subset$Phenotypic_Class <- factor(\n",
    "            coef_subset$Phenotypic_Class,\n",
    "            levels=unique(coef_subset$Phenotypic_Class)\n",
    "        )\n",
    "        coef_subset_mat <- as.matrix(coef_subset[, -1])\n",
    "        rownames(coef_subset_mat) <- coef_subset$Phenotypic_Class\n",
    "\n",
    "        # Process metadata for inclusion in heatmap annotation\n",
    "        metadata_subset <- dplyr::as_tibble(\n",
    "            colnames(coef_subset_mat),\n",
    "            .name_repair = function(x) \"feature\"\n",
    "            ) %>%\n",
    "            tidyr::separate(\n",
    "                feature,\n",
    "                into = first_component_split,\n",
    "                sep = \"__\",\n",
    "                remove = FALSE\n",
    "            ) %>%\n",
    "            tidyr::separate(\n",
    "                \"feature_name\",\n",
    "                into = metadata_annotation_options[[feature_type_value]],\n",
    "                sep = \"_\",\n",
    "                remove = FALSE\n",
    "            )\n",
    "        metadata_subset <- as.matrix(metadata_subset)\n",
    "        rownames(metadata_subset) <- colnames(coef_subset_mat)\n",
    "        \n",
    "        # Store in list\n",
    "        feature_matrix_list[[list_name_index]] <- list()\n",
    "        feature_matrix_list[[list_name_index]][[\"coef_matrix\"]] <- coef_subset_mat\n",
    "        feature_matrix_list[[list_name_index]][[\"min_max_for_plotting\"]] <- c(min(coef_subset_mat), max(coef_subset_mat))\n",
    "        feature_matrix_list[[list_name_index]][[\"metadata_annotation\"]] <- metadata_subset\n",
    "    }\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "8ad22f62-5f64-4d55-85a8-4a60b7a45f4d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] \"Generating heatmap for: CP__False\"\n",
      "[1] \"Generating heatmap for: CP_and_DP__False\"\n",
      "[1] \"Generating heatmap for: DP__False\"\n",
      "[1] \"Generating heatmap for: CP__True\"\n",
      "[1] \"Generating heatmap for: CP_and_DP__True\"\n",
      "[1] \"Generating heatmap for: DP__True\"\n"
     ]
    }
   ],
   "source": [
    "model_heatmap_file_names <- list()\n",
    "for (model in names(feature_matrix_list)) {\n",
    "    print(paste(\"Generating heatmap for:\", model))\n",
    "\n",
    "    # Create output file\n",
    "    output_file <- file.path(output_dir, paste0(\"heatmap_\", model, \".png\"))\n",
    "\n",
    "    # Store model file in a list for downstream loading\n",
    "    model_heatmap_file_names[model] <- output_file\n",
    "    \n",
    "    # Create components for plotheatmap_gg_list[[\"CP__False\"]]ting and subsetting\n",
    "    model_split_details <- unlist(stringr::str_split(model, \"__\"))\n",
    "    feature_space <- model_split_details[1]\n",
    "    shuffled_or_not <- model_split_details[2]\n",
    "\n",
    "    if (shuffled_or_not == \"False\") {\n",
    "        column_title = \"Real data (final model)\"\n",
    "    } else {\n",
    "        column_title = \"Shuffled data\"\n",
    "    }\n",
    "\n",
    "    # Generate heatmaps depending on the feature space\n",
    "    if (feature_space == \"CP\") {\n",
    "        column_title <- paste(\"CellProfiler features\", column_title)\n",
    "        heatmap_color_range <- c(\n",
    "            feature_matrix_list[[\"CP__False\"]][[\"min_max_for_plotting\"]][1], 0,\n",
    "            feature_matrix_list[[\"CP__False\"]][[\"min_max_for_plotting\"]][2]\n",
    "        )\n",
    "        \n",
    "        coef_heatmap <- Heatmap(\n",
    "            feature_matrix_list[[model]][[\"coef_matrix\"]],\n",
    "            col = circlize::colorRamp2(heatmap_color_range, c(\"blue\", \"white\", \"red\")),\n",
    "            top_annotation = HeatmapAnnotation(\n",
    "                df = as.data.frame(feature_matrix_list[[model]][[\"metadata_annotation\"]]) %>%\n",
    "                    dplyr::select(feature_group),\n",
    "                col = list(feature_group = cp_feature_group_colors),\n",
    "                annotation_legend_param = list(feature_group = list(title = \"CP feature\\ngroup\"))\n",
    "            ),\n",
    "            column_split = as.data.frame(feature_matrix_list[[model]][[\"metadata_annotation\"]])$feature_group,\n",
    "            column_title = column_title,\n",
    "            name = \"ML Coefficient\",\n",
    "            show_column_names = FALSE\n",
    "        )\n",
    "\n",
    "    } else if (feature_space == \"CP_and_DP\") {\n",
    "        column_title <- paste(\"CP and DP features\", column_title)\n",
    "        heatmap_color_range <- c(\n",
    "            feature_matrix_list[[\"CP_and_DP__False\"]][[\"min_max_for_plotting\"]][1], 0,\n",
    "            feature_matrix_list[[\"CP_and_DP__False\"]][[\"min_max_for_plotting\"]][2]\n",
    "        )\n",
    "        \n",
    "        coef_heatmap <- Heatmap(\n",
    "            feature_matrix_list[[model]][[\"coef_matrix\"]],\n",
    "            col = circlize::colorRamp2(heatmap_color_range, c(\"blue\", \"white\", \"red\")),\n",
    "            top_annotation = HeatmapAnnotation(\n",
    "                df = as.data.frame(feature_matrix_list[[model]][[\"metadata_annotation\"]]) %>%\n",
    "                    dplyr::select(feature_space),\n",
    "                col = list(feature_space = feature_space_colors),\n",
    "                annotation_legend_param = list(feature_space = list(title = \"Feature\\nspace\"))\n",
    "            ),\n",
    "            column_title = column_title,\n",
    "            name = \"ML Coefficient\",\n",
    "            show_column_names = FALSE\n",
    "        )\n",
    "    } else {\n",
    "        column_title <- paste(\"DeepProfiler features\", column_title)\n",
    "        heatmap_color_range <- c(\n",
    "            feature_matrix_list[[\"DP__False\"]][[\"min_max_for_plotting\"]][1], 0,\n",
    "            feature_matrix_list[[\"DP__False\"]][[\"min_max_for_plotting\"]][2]\n",
    "        )\n",
    "        \n",
    "        coef_heatmap <- Heatmap(\n",
    "            feature_matrix_list[[model]][[\"coef_matrix\"]],\n",
    "            col = circlize::colorRamp2(heatmap_color_range, c(\"blue\", \"white\", \"red\")),\n",
    "            top_annotation = HeatmapAnnotation(\n",
    "                df = as.data.frame(feature_matrix_list[[model]][[\"metadata_annotation\"]]) %>%\n",
    "                    dplyr::select(feature_space),\n",
    "                col = list(feature_space = feature_space_colors),\n",
    "                annotation_legend_param = list(feature_space = list(title = \"Feature\\nspace\"))\n",
    "            ),\n",
    "            column_split = as.data.frame(feature_matrix_list[[model]][[\"metadata_annotation\"]])$feature_space,\n",
    "            column_title = column_title,\n",
    "            name = \"ML Coefficient\",\n",
    "            show_column_names = FALSE\n",
    "        )\n",
    "    }\n",
    "\n",
    "    png(output_file, width = 570, height = 275)\n",
    "    draw(coef_heatmap, merge_legend = TRUE)\n",
    "    dev.off()\n",
    "}"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "4.2.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}