3.processing_features/2.pycytominer_singlecell_pipelines.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "jukit_cell_id": "kVVzyv0WPc"
   },
   "source": [
    "## Perform single-cell pycytominer pipelines\n",
    "\n",
    "Following single-cell curation with cytotable, we create single-cell profiles by applying the following steps:\n",
    "\n",
    "1. annotation\n",
    "2. normalization\n",
    "3. feature_selection\n",
    "\n",
    "Additionally, we create bulk profiles following feature selection.\n",
    "We call this \"Cameron's Method\".\n",
    "\n",
    "4. Aggregate (to form bulk, after single-cell processing)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-03-21T18:57:16.153188Z",
     "iopub.status.busy": "2024-03-21T18:57:16.153072Z",
     "iopub.status.idle": "2024-03-21T18:57:17.159989Z",
     "shell.execute_reply": "2024-03-21T18:57:17.159635Z"
    },
    "jukit_cell_id": "wow8TtvKeA"
   },
   "outputs": [],
   "source": [
    "import pathlib\n",
    "import yaml\n",
    "import pprint\n",
    "\n",
    "import pandas as pd\n",
    "\n",
    "from pycytominer import aggregate, annotate, normalize, feature_select\n",
    "from pycytominer.cyto_utils import load_profiles, output, infer_cp_features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-03-21T18:57:17.162051Z",
     "iopub.status.busy": "2024-03-21T18:57:17.161892Z",
     "iopub.status.idle": "2024-03-21T18:57:17.165965Z",
     "shell.execute_reply": "2024-03-21T18:57:17.165686Z"
    },
    "jukit_cell_id": "UUkCcbzcrv"
   },
   "outputs": [],
   "source": [
    "# Set constants\n",
    "feature_select_ops = [\n",
    "    \"variance_threshold\",\n",
    "    \"correlation_threshold\",\n",
    "    \"blocklist\",\n",
    "    \"drop_na_columns\",\n",
    "]\n",
    "\n",
    "# Columns to remove prior to single-cell aggregation via cameron's method\n",
    "cameron_unwanted_aggregate_cols = {\"Object\", \"Parent\", \"Site\", \"Image\", \"Location\"}\n",
    "\n",
    "# Set paths\n",
    "output_dir = pathlib.Path(\"data/single_cell_profiles\")\n",
    "output_dir.mkdir(exist_ok=True)\n",
    "metadata_dir = pathlib.Path(\"../0.download_data/metadata/\")\n",
    "\n",
    "# load in plate information\n",
    "dictionary_path = pathlib.Path(\"./plate_info_dictionary.yaml\")\n",
    "with open(dictionary_path) as file:\n",
    "    plate_info_dictionary = yaml.load(file, Loader=yaml.FullLoader)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-03-21T18:57:17.167330Z",
     "iopub.status.busy": "2024-03-21T18:57:17.167224Z",
     "iopub.status.idle": "2024-03-21T18:57:17.170649Z",
     "shell.execute_reply": "2024-03-21T18:57:17.170423Z"
    },
    "jukit_cell_id": "NclnZZuqrJ"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{   'Plate_3': {   'dest_path': 'data/converted_data/Plate_3.parquet',\n",
      "                   'platemap_path': '/home/jenna/nf1_cellpainting_data/0.download_data/metadata/platemap_NF1_plate3.csv',\n",
      "                   'source_path': '/home/jenna/nf1_cellpainting_data/2.cellprofiler_analysis/analysis_output/Plate_3/Plate_3_nf1_analysis.sqlite'},\n",
      "    'Plate_3_prime': {   'dest_path': 'data/converted_data/Plate_3_prime.parquet',\n",
      "                         'platemap_path': '/home/jenna/nf1_cellpainting_data/0.download_data/metadata/platemap_NF1_plate3.csv',\n",
      "                         'source_path': '/home/jenna/nf1_cellpainting_data/2.cellprofiler_analysis/analysis_output/Plate_3_prime/Plate_3_prime_nf1_analysis.sqlite'}}\n"
     ]
    }
   ],
   "source": [
    "# add path to platemaps for each plate\n",
    "for plate in plate_info_dictionary.keys():\n",
    "    # since Plate_3_prime has the same platemap as Plate_3,\n",
    "    # we need an else statement so that we make sure it adds the\n",
    "    # path that was given to Plate_3\n",
    "    if plate != \"Plate_3_prime\":\n",
    "        # match the naming format of the plates to the platemap file\n",
    "        plate_info_dictionary[plate][\"platemap_path\"] = str(\n",
    "            pathlib.Path(\n",
    "                list(\n",
    "                    metadata_dir.rglob(\n",
    "                        f\"platemap_NF1_{plate.replace('_', '').lower()}.csv\"\n",
    "                    )\n",
    "                )[0]\n",
    "            ).resolve(strict=True)\n",
    "        )\n",
    "    else:\n",
    "        plate_info_dictionary[\"Plate_3_prime\"][\"platemap_path\"] = plate_info_dictionary[\n",
    "            \"Plate_3\"\n",
    "        ][\"platemap_path\"]\n",
    "\n",
    "# view the dictionary to assess that all info is added correctly\n",
    "pprint.pprint(plate_info_dictionary, indent=4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "jukit_cell_id": "EZ5YZymnlX"
   },
   "source": [
    "## Perform single-cell pycytominer pipeline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-03-21T18:57:17.190993Z",
     "iopub.status.busy": "2024-03-21T18:57:17.190677Z",
     "iopub.status.idle": "2024-03-21T19:00:43.498737Z",
     "shell.execute_reply": "2024-03-21T19:00:43.498183Z"
    },
    "jukit_cell_id": "meValC3kNF"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Now performing single-cell pycytominer pipeline for Plate_3\n",
      "HET cells have been removed from Plate_3\n",
      "Performing normalization for Plate_3 using samples parameter: all\n",
      "(48, 1161)\n",
      "Now performing single-cell pycytominer pipeline for Plate_3_prime\n",
      "HET cells have been removed from Plate_3_prime\n",
      "Performing normalization for Plate_3_prime using samples parameter: all\n",
      "(48, 1136)\n"
     ]
    }
   ],
   "source": [
    "for plate, info in plate_info_dictionary.items():\n",
    "    print(f\"Now performing single-cell pycytominer pipeline for {plate}\")\n",
    "    output_annotated_file = str(\n",
    "        pathlib.Path(f\"{output_dir}/{plate}_sc_annotated.parquet\")\n",
    "    )\n",
    "    output_normalized_file = str(\n",
    "        pathlib.Path(f\"{output_dir}/{plate}_sc_normalized.parquet\")\n",
    "    )\n",
    "    output_feature_select_file = str(\n",
    "        pathlib.Path(f\"{output_dir}/{plate}_sc_feature_selected.parquet\")\n",
    "    )\n",
    "    output_aggregated_file = str(\n",
    "        pathlib.Path(f\"{output_dir}/{plate}_bulk_camerons_method.parquet\")\n",
    "    )\n",
    "\n",
    "    # Load single-cell profiles\n",
    "    single_cell_df = pd.read_parquet(info[\"dest_path\"])\n",
    "\n",
    "    # Load plate map\n",
    "    platemap_df = pd.read_csv(info[\"platemap_path\"])\n",
    "\n",
    "    # Step 1: Annotation\n",
    "    # add metadata from plate map file to extracted single cell features\n",
    "    annotated_df = annotate(\n",
    "        profiles=single_cell_df,\n",
    "        platemap=platemap_df,\n",
    "        join_on=[\"Metadata_well_position\", \"Image_Metadata_Well\"],\n",
    "    )\n",
    "\n",
    "    # rename site column to avoid any issues with identifying the column as metadata over feature\n",
    "    annotated_df = annotated_df.rename(columns={\"Image_Metadata_Site\": \"Metadata_Site\"})\n",
    "\n",
    "    # move metadata well, single cell count, and site to the front of the df (for easy visualization in python)\n",
    "    well_column = annotated_df.pop(\"Metadata_Well\")\n",
    "    singlecell_column = annotated_df.pop(\"Metadata_number_of_singlecells\")\n",
    "    site_column = annotated_df.pop(\"Metadata_Site\")\n",
    "\n",
    "    # insert the columns in specific parts of the data frame\n",
    "    annotated_df.insert(2, \"Metadata_Well\", well_column)\n",
    "    annotated_df.insert(3, \"Metadata_Site\", site_column)\n",
    "    annotated_df.insert(4, \"Metadata_number_of_singlecells\", singlecell_column)\n",
    "\n",
    "    # for plates 5, 3, and 3 prime, remove any rows with HET due to contamination/not using in model\n",
    "    if plate in [\"Plate_5\", \"Plate_3\", \"Plate_3_prime\"]:\n",
    "        # Filter single-cell profiles, removing HET genotype\n",
    "        annotated_df = annotated_df[annotated_df[\"Metadata_genotype\"] != \"HET\"].reset_index(drop=True)\n",
    "        print(\"HET cells have been removed from\", plate)\n",
    "\n",
    "    # use output to save annotated df as you can not use the output parameter in the annotate which will not return data frame\n",
    "    output(\n",
    "        df=annotated_df,\n",
    "        output_filename=output_annotated_file,\n",
    "        output_type=\"parquet\",\n",
    "    )\n",
    "\n",
    "    # set default for samples to use in normalization\n",
    "    samples = \"all\"\n",
    "\n",
    "    # Only for Plate 4, we want to normalize to no siRNA treatment Null and WT cells (controls)\n",
    "    if plate == \"Plate_4\":\n",
    "        samples = \"Metadata_Concentration == 0.0 and (Metadata_genotype == 'Null' or Metadata_genotype == 'WT')\"\n",
    "\n",
    "\n",
    "    print(f\"Performing normalization for {plate} using samples parameter: {samples}\")\n",
    "\n",
    "    # Step 2: Normalization\n",
    "    normalized_df = normalize(\n",
    "        profiles=output_annotated_file,\n",
    "        method=\"standardize\",\n",
    "        output_file=output_normalized_file,\n",
    "        output_type=\"parquet\",\n",
    "        samples=samples,\n",
    "    )\n",
    "\n",
    "    # Step 3: Feature selection\n",
    "    feature_select(\n",
    "        output_normalized_file,\n",
    "        operation=feature_select_ops,\n",
    "        na_cutoff=0,\n",
    "        output_file=output_feature_select_file,\n",
    "        output_type=\"parquet\",\n",
    "    )\n",
    "\n",
    "    # Step 4: Cameron's method of aggregation\n",
    "    feature_select_df = load_profiles(output_feature_select_file)\n",
    "    # Specify metadata columns in aggregation step to ensure they are retained for downstream analysis\n",
    "    metadata_cols = infer_cp_features(feature_select_df, metadata=True)\n",
    "    metadata_cols = [\n",
    "        x\n",
    "        for x in metadata_cols\n",
    "        if all(col not in x for col in cameron_unwanted_aggregate_cols)\n",
    "    ]\n",
    "\n",
    "    aggregate_df = aggregate(\n",
    "        population_df=feature_select_df,\n",
    "        operation=\"median\",\n",
    "        strata=metadata_cols,\n",
    "        output_file=output_aggregated_file,\n",
    "        output_type=\"parquet\",\n",
    "    )\n",
    "\n",
    "    print(aggregate_df.shape)"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "python",
   "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.8.18"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}