notebooks/mitocheck-map-analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# MAP Analysis with MitoCheck Single Cells \n",
    "\n",
    "\n",
    "In this notebook, our goal is to apply the Mean Average Precision (MAP) metric, developed in the [copairs](https://github.com/cytomining/copairs) analysis package.\n",
    "We apply this metric to the MitoCheck single-cell dataset, to see the effects of genetic pertubations based on their phenotype. \n",
    "\n",
    "\n",
    "Some links to look at:\n",
    "- copairs [repo](https://github.com/cytomining/copairs)\n",
    "- MitoCheck github [repo](https://github.com/WayScience/mitocheck_data)\n",
    "- MitoCheck zenodo [repo](https://zenodo.org/records/7967386)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "import logging\n",
    "import pprint\n",
    "import pathlib\n",
    "\n",
    "from copairs.map import run_pipeline\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "\n",
    "# imports src\n",
    "sys.path.append(\"../\")\n",
    "from src import utils\n",
    "\n",
    "# setting up logger\n",
    "# setting up logger\n",
    "logging.basicConfig(filename=\"map_analysis.log\", \n",
    "                    level=logging.DEBUG,\n",
    "                    format='%(levelname)s:%(asctime)s:%(name)s:%(message)s')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Loading Downloaded Data\n",
    "\n",
    "In this section, we load the MitoCheck single-cell datasets, including the training, positive controls, and negative controls. \n",
    "For detailed information about the dataset, please refer to the MitoCheck report mentioned above.\n",
    "\n",
    "After downloading the data, we perform formatting by dividing it into two sections. \n",
    "The first section comprises the metadata of each individual cell, while the second section presents all quantified features in a numpy array format.\n",
    "\n",
    "This formatting is designed to easily integrate with the copairs `run_pipeline()` function, allowing for easy execution of the analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# parameters\n",
    "training_singlecell_data = pathlib.Path(\"../data/raw/training_data.csv.gz\").resolve(strict=True)\n",
    "pos_control_data = pathlib.Path(\"../data/raw/normalized_data/positive_control_data.csv.gz\").resolve(strict=True)\n",
    "neg_control_data = pathlib.Path(\"../data/raw/normalized_data/negative_control_data.csv.gz\").resolve(strict=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# # loading in the data into dataframe (~10min loading)\n",
    "# training_sc_data = pd.read_csv(training_singlecell_data).drop(\"Unnamed: 0\", axis=1)\n",
    "# pos_control_sc_data = pd.read_csv(pos_control_data)\n",
    "# neg_control_sc_data = pd.read_csv(neg_control_data)\n",
    "\n",
    "# # adding the Mitocheck_Phenotypic_Class into the controls  and labels\n",
    "# pos_control_sc_data.insert(0, \"Mitocheck_Phenotypic_Class\", \"pos_control\")\n",
    "# neg_control_sc_data.insert(0, \"Mitocheck_Phenotypic_Class\", \"neg_control\")\n",
    "\n",
    "# # droping column from trainign data since it does not exist in the controls \n",
    "# training_sc_data = training_sc_data.drop(\"Metadata_Object_Outline\", axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO: delete this later for \n",
    "training_sc_data = pd.read_parquet(\"../data/processed/training_sc_data.parquet\")\n",
    "pos_control_sc_data = pd.read_parquet(\"../data/processed/pos_control_sc_data.parquet\").sample(frac=0.01, random_state=42)\n",
    "neg_control_sc_data = pd.read_parquet(\"../data/processed/neg_control_sc_data.parquet\").sample(frac=0.01, random_state=42)\n",
    "\n",
    "# adding the Mitocheck_Phenotypic_Class into the controls  and labels\n",
    "neg_control_sc_data.insert(0, \"Mitocheck_Phenotypic_Class\", \"neg_control\")\n",
    "\n",
    "# adding control labels into the dataset\n",
    "training_sc_data.insert(1, \"Metadata_is_control\", 0)\n",
    "neg_control_sc_data.insert(1, \"Metadata_is_control\", 1)\n",
    "\n",
    "# droping column from trainign data since it does not exist in the controls \n",
    "training_sc_data = training_sc_data.drop(\"Metadata_Object_Outline\", axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "ename": "KeyError",
     "evalue": "\"['Metadata_is_control == 0'] not in index\"",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mKeyError\u001b[0m                                  Traceback (most recent call last)",
      "\u001b[1;32m/home/erikserrano/Development/Mitocheck-MAP-analysis/notebooks/mitocheck-map-analysis.ipynb Cell 7\u001b[0m line \u001b[0;36m3\n\u001b[1;32m     <a href='vscode-notebook-cell:/home/erikserrano/Development/Mitocheck-MAP-analysis/notebooks/mitocheck-map-analysis.ipynb#W6sZmlsZQ%3D%3D?line=32'>33</a>\u001b[0m \u001b[39m# execute pipeline on negative control with trianing dataset with cp features\u001b[39;00m\n\u001b[1;32m     <a href='vscode-notebook-cell:/home/erikserrano/Development/Mitocheck-MAP-analysis/notebooks/mitocheck-map-analysis.ipynb#W6sZmlsZQ%3D%3D?line=33'>34</a>\u001b[0m logging\u001b[39m.\u001b[39minfo(\u001b[39mf\u001b[39m\u001b[39m\"\u001b[39m\u001b[39mRunning pipeline on CP features using \u001b[39m\u001b[39m{\u001b[39;00mphenotype\u001b[39m}\u001b[39;00m\u001b[39m phenotype\u001b[39m\u001b[39m\"\u001b[39m)\n\u001b[0;32m---> <a href='vscode-notebook-cell:/home/erikserrano/Development/Mitocheck-MAP-analysis/notebooks/mitocheck-map-analysis.ipynb#W6sZmlsZQ%3D%3D?line=34'>35</a>\u001b[0m cp_negative_training_result \u001b[39m=\u001b[39m run_pipeline(meta\u001b[39m=\u001b[39;49mnegative_training_cp_meta,\n\u001b[1;32m     <a href='vscode-notebook-cell:/home/erikserrano/Development/Mitocheck-MAP-analysis/notebooks/mitocheck-map-analysis.ipynb#W6sZmlsZQ%3D%3D?line=35'>36</a>\u001b[0m                                         feats\u001b[39m=\u001b[39;49mnegative_training_cp_feats,\n\u001b[1;32m     <a href='vscode-notebook-cell:/home/erikserrano/Development/Mitocheck-MAP-analysis/notebooks/mitocheck-map-analysis.ipynb#W6sZmlsZQ%3D%3D?line=36'>37</a>\u001b[0m                                         pos_sameby\u001b[39m=\u001b[39;49mpos_sameby,\n\u001b[1;32m     <a href='vscode-notebook-cell:/home/erikserrano/Development/Mitocheck-MAP-analysis/notebooks/mitocheck-map-analysis.ipynb#W6sZmlsZQ%3D%3D?line=37'>38</a>\u001b[0m                                         pos_diffby\u001b[39m=\u001b[39;49mpos_diffby,\n\u001b[1;32m     <a href='vscode-notebook-cell:/home/erikserrano/Development/Mitocheck-MAP-analysis/notebooks/mitocheck-map-analysis.ipynb#W6sZmlsZQ%3D%3D?line=38'>39</a>\u001b[0m                                         neg_sameby\u001b[39m=\u001b[39;49mneg_sameby,\n\u001b[1;32m     <a href='vscode-notebook-cell:/home/erikserrano/Development/Mitocheck-MAP-analysis/notebooks/mitocheck-map-analysis.ipynb#W6sZmlsZQ%3D%3D?line=39'>40</a>\u001b[0m                                         neg_diffby\u001b[39m=\u001b[39;49mneg_diffby,\n\u001b[1;32m     <a href='vscode-notebook-cell:/home/erikserrano/Development/Mitocheck-MAP-analysis/notebooks/mitocheck-map-analysis.ipynb#W6sZmlsZQ%3D%3D?line=40'>41</a>\u001b[0m                                         batch_size\u001b[39m=\u001b[39;49mbatch_size,\n\u001b[1;32m     <a href='vscode-notebook-cell:/home/erikserrano/Development/Mitocheck-MAP-analysis/notebooks/mitocheck-map-analysis.ipynb#W6sZmlsZQ%3D%3D?line=41'>42</a>\u001b[0m                                         null_size\u001b[39m=\u001b[39;49mnull_size)\n\u001b[1;32m     <a href='vscode-notebook-cell:/home/erikserrano/Development/Mitocheck-MAP-analysis/notebooks/mitocheck-map-analysis.ipynb#W6sZmlsZQ%3D%3D?line=42'>43</a>\u001b[0m map_results_neg_cp\u001b[39m.\u001b[39mappend(cp_negative_training_result)                                       \n\u001b[1;32m     <a href='vscode-notebook-cell:/home/erikserrano/Development/Mitocheck-MAP-analysis/notebooks/mitocheck-map-analysis.ipynb#W6sZmlsZQ%3D%3D?line=44'>45</a>\u001b[0m \u001b[39m# execute pipeline on negative control with trianing dataset with dp features\u001b[39;00m\n",
      "File \u001b[0;32m~/Programs/miniconda3/envs/map/lib/python3.12/site-packages/copairs/map.py:115\u001b[0m, in \u001b[0;36mrun_pipeline\u001b[0;34m(meta, feats, pos_sameby, pos_diffby, neg_sameby, neg_diffby, null_size, batch_size, seed)\u001b[0m\n\u001b[1;32m    105\u001b[0m \u001b[39mdef\u001b[39;00m \u001b[39mrun_pipeline\u001b[39m(meta,\n\u001b[1;32m    106\u001b[0m                  feats,\n\u001b[1;32m    107\u001b[0m                  pos_sameby,\n\u001b[0;32m   (...)\u001b[0m\n\u001b[1;32m    112\u001b[0m                  batch_size\u001b[39m=\u001b[39m\u001b[39m20000\u001b[39m,\n\u001b[1;32m    113\u001b[0m                  seed\u001b[39m=\u001b[39m\u001b[39m0\u001b[39m) \u001b[39m-\u001b[39m\u001b[39m>\u001b[39m pd\u001b[39m.\u001b[39mDataFrame:\n\u001b[1;32m    114\u001b[0m     columns \u001b[39m=\u001b[39m flatten_str_list(pos_sameby, pos_diffby, neg_sameby, neg_diffby)\n\u001b[0;32m--> 115\u001b[0m     validate_pipeline_input(meta, feats, columns)\n\u001b[1;32m    117\u001b[0m     \u001b[39m# Critical!, otherwise the indexing wont work\u001b[39;00m\n\u001b[1;32m    118\u001b[0m     meta \u001b[39m=\u001b[39m meta\u001b[39m.\u001b[39mreset_index(drop\u001b[39m=\u001b[39m\u001b[39mTrue\u001b[39;00m)\u001b[39m.\u001b[39mcopy()\n",
      "File \u001b[0;32m~/Programs/miniconda3/envs/map/lib/python3.12/site-packages/copairs/map.py:99\u001b[0m, in \u001b[0;36mvalidate_pipeline_input\u001b[0;34m(meta, feats, columns)\u001b[0m\n\u001b[1;32m     98\u001b[0m \u001b[39mdef\u001b[39;00m \u001b[39mvalidate_pipeline_input\u001b[39m(meta, feats, columns):\n\u001b[0;32m---> 99\u001b[0m     \u001b[39mif\u001b[39;00m meta[columns]\u001b[39m.\u001b[39misna()\u001b[39m.\u001b[39many(axis\u001b[39m=\u001b[39m\u001b[39mNone\u001b[39;00m):\n\u001b[1;32m    100\u001b[0m         \u001b[39mraise\u001b[39;00m \u001b[39mValueError\u001b[39;00m(\u001b[39m'\u001b[39m\u001b[39mmetadata columns should not have null values.\u001b[39m\u001b[39m'\u001b[39m)\n\u001b[1;32m    101\u001b[0m     \u001b[39mif\u001b[39;00m \u001b[39mlen\u001b[39m(meta) \u001b[39m!=\u001b[39m \u001b[39mlen\u001b[39m(feats):\n",
      "File \u001b[0;32m~/Programs/miniconda3/envs/map/lib/python3.12/site-packages/pandas/core/frame.py:3899\u001b[0m, in \u001b[0;36mDataFrame.__getitem__\u001b[0;34m(self, key)\u001b[0m\n\u001b[1;32m   3897\u001b[0m     \u001b[39mif\u001b[39;00m is_iterator(key):\n\u001b[1;32m   3898\u001b[0m         key \u001b[39m=\u001b[39m \u001b[39mlist\u001b[39m(key)\n\u001b[0;32m-> 3899\u001b[0m     indexer \u001b[39m=\u001b[39m \u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49mcolumns\u001b[39m.\u001b[39;49m_get_indexer_strict(key, \u001b[39m\"\u001b[39;49m\u001b[39mcolumns\u001b[39;49m\u001b[39m\"\u001b[39;49m)[\u001b[39m1\u001b[39m]\n\u001b[1;32m   3901\u001b[0m \u001b[39m# take() does not accept boolean indexers\u001b[39;00m\n\u001b[1;32m   3902\u001b[0m \u001b[39mif\u001b[39;00m \u001b[39mgetattr\u001b[39m(indexer, \u001b[39m\"\u001b[39m\u001b[39mdtype\u001b[39m\u001b[39m\"\u001b[39m, \u001b[39mNone\u001b[39;00m) \u001b[39m==\u001b[39m \u001b[39mbool\u001b[39m:\n",
      "File \u001b[0;32m~/Programs/miniconda3/envs/map/lib/python3.12/site-packages/pandas/core/indexes/base.py:6114\u001b[0m, in \u001b[0;36mIndex._get_indexer_strict\u001b[0;34m(self, key, axis_name)\u001b[0m\n\u001b[1;32m   6111\u001b[0m \u001b[39melse\u001b[39;00m:\n\u001b[1;32m   6112\u001b[0m     keyarr, indexer, new_indexer \u001b[39m=\u001b[39m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39m_reindex_non_unique(keyarr)\n\u001b[0;32m-> 6114\u001b[0m \u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49m_raise_if_missing(keyarr, indexer, axis_name)\n\u001b[1;32m   6116\u001b[0m keyarr \u001b[39m=\u001b[39m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mtake(indexer)\n\u001b[1;32m   6117\u001b[0m \u001b[39mif\u001b[39;00m \u001b[39misinstance\u001b[39m(key, Index):\n\u001b[1;32m   6118\u001b[0m     \u001b[39m# GH 42790 - Preserve name from an Index\u001b[39;00m\n",
      "File \u001b[0;32m~/Programs/miniconda3/envs/map/lib/python3.12/site-packages/pandas/core/indexes/base.py:6178\u001b[0m, in \u001b[0;36mIndex._raise_if_missing\u001b[0;34m(self, key, indexer, axis_name)\u001b[0m\n\u001b[1;32m   6175\u001b[0m     \u001b[39mraise\u001b[39;00m \u001b[39mKeyError\u001b[39;00m(\u001b[39mf\u001b[39m\u001b[39m\"\u001b[39m\u001b[39mNone of [\u001b[39m\u001b[39m{\u001b[39;00mkey\u001b[39m}\u001b[39;00m\u001b[39m] are in the [\u001b[39m\u001b[39m{\u001b[39;00maxis_name\u001b[39m}\u001b[39;00m\u001b[39m]\u001b[39m\u001b[39m\"\u001b[39m)\n\u001b[1;32m   6177\u001b[0m not_found \u001b[39m=\u001b[39m \u001b[39mlist\u001b[39m(ensure_index(key)[missing_mask\u001b[39m.\u001b[39mnonzero()[\u001b[39m0\u001b[39m]]\u001b[39m.\u001b[39munique())\n\u001b[0;32m-> 6178\u001b[0m \u001b[39mraise\u001b[39;00m \u001b[39mKeyError\u001b[39;00m(\u001b[39mf\u001b[39m\u001b[39m\"\u001b[39m\u001b[39m{\u001b[39;00mnot_found\u001b[39m}\u001b[39;00m\u001b[39m not in index\u001b[39m\u001b[39m\"\u001b[39m)\n",
      "\u001b[0;31mKeyError\u001b[0m: \"['Metadata_is_control == 0'] not in index\""
     ]
    }
   ],
   "source": [
    "# parameters for pipeline\n",
    "pos_sameby = [\"Mitocheck_Phenotypic_Class\", \"Metadata_is_control == 0\"]\n",
    "pos_diffby = []\n",
    "neg_sameby = [\"Metadata_Plate\"]\n",
    "neg_diffby = [\"Metadata_is_control\"]\n",
    "null_size = 3000\n",
    "batch_size = 3000\n",
    "\n",
    "# storing all map results based on postiive and negative controls and feature types\n",
    "map_results_neg_cp = []\n",
    "map_results_neg_dp = []\n",
    "map_results_neg_cp_dp = []\n",
    "\n",
    "# running process\n",
    "# for loop selects one single phenotype\n",
    "# then splits the data into metadata and raw feature values \n",
    "# two different groups that contains 3 splits caused by the types of features\n",
    "# applie the copairs pipeline\n",
    "for phenotype in list(training_sc_data[\"Mitocheck_Phenotypic_Class\"].unique()):\n",
    "\n",
    "    # select training dataset based on phenotype\n",
    "    selected_training = training_sc_data.loc[training_sc_data[\"Mitocheck_Phenotypic_Class\"] == phenotype]\n",
    "\n",
    "    # concatenate to positive and negative control \n",
    "    training_w_neg = pd.concat([selected_training, neg_control_sc_data])\n",
    "    \n",
    "    # spliting metadata and raw feature values\n",
    "    negative_training_cp_meta, negative_training_cp_feats = utils.split_data(training_w_neg, dataset=\"CP\")\n",
    "    negative_training_dp_meta, negative_training_dp_feats = utils.split_data(training_w_neg, dataset=\"DP\")\n",
    "    negative_training_cp_dp_meta, negative_training_cp_dp_feats = utils.split_data(training_w_neg, dataset=\"CP_and_DP\")\n",
    "\n",
    "\n",
    "    # execute pipeline on negative control with trianing dataset with cp features\n",
    "    logging.info(f\"Running pipeline on CP features using {phenotype} phenotype\")\n",
    "    cp_negative_training_result = run_pipeline(meta=negative_training_cp_meta,\n",
    "                                            feats=negative_training_cp_feats,\n",
    "                                            pos_sameby=pos_sameby,\n",
    "                                            pos_diffby=pos_diffby,\n",
    "                                            neg_sameby=neg_sameby,\n",
    "                                            neg_diffby=neg_diffby,\n",
    "                                            batch_size=batch_size,\n",
    "                                            null_size=null_size)\n",
    "    map_results_neg_cp.append(cp_negative_training_result)                                       \n",
    "\n",
    "    # execute pipeline on negative control with trianing dataset with dp features\n",
    "    logging.info(f\"Running pipeline on DP features using {phenotype} phenotype\")\n",
    "    dp_negative_training_result = run_pipeline(meta=negative_training_dp_meta,\n",
    "                                            feats=negative_training_dp_feats,\n",
    "                                            pos_sameby=pos_sameby,\n",
    "                                            pos_diffby=pos_diffby,\n",
    "                                            neg_sameby=neg_sameby,\n",
    "                                            neg_diffby=neg_diffby,\n",
    "                                            batch_size=batch_size,\n",
    "                                            null_size=null_size)\n",
    "    map_results_neg_dp.append(dp_negative_training_result)                                       \n",
    "\n",
    "    # execute pipeline on negative control with trianing dataset with cp_dp features\n",
    "    logging.info(f\"Running pipeline on CP and DP features using {phenotype} phenotype\")\n",
    "    cp_dp_negative_training_result = run_pipeline(meta=negative_training_cp_dp_meta,\n",
    "                                            feats=negative_training_cp_dp_feats,\n",
    "                                            pos_sameby=pos_sameby,\n",
    "                                            pos_diffby=pos_diffby,\n",
    "                                            neg_sameby=neg_sameby,\n",
    "                                            neg_diffby=neg_diffby,\n",
    "                                            batch_size=batch_size,\n",
    "                                            null_size=null_size)\n",
    "    map_results_neg_cp_dp.append(cp_dp_negative_training_result)                                       "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "mitocheck-map-analysis",
   "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.12.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}