3.processing_features/1.pycytominer_bulk_pipelines.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "91a846b5-8ead-4ad8-b291-00c09a09ab0a",
   "metadata": {},
   "source": [
    "## Perform traditional bulk pycytominer pipeline\n",
    "\n",
    "Following single-cell curation with cytotable, we create bulk profiles by applying the following steps:\n",
    "\n",
    "1. aggregation\n",
    "2. annotation\n",
    "3. normalization\n",
    "4. feature_selection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "4b02e669-6b8b-49d7-a0a0-678eff76a854",
   "metadata": {},
   "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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "609183fc-b780-4731-9a63-cbfb448d4f66",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set constants\n",
    "feature_select_ops = [\n",
    "    \"variance_threshold\",\n",
    "    \"correlation_threshold\",\n",
    "    \"blocklist\",\n",
    "    \"drop_na_columns\",\n",
    "]\n",
    "\n",
    "# Set paths\n",
    "output_dir = pathlib.Path(\"data/bulk_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,
   "id": "88811904-baa5-47a4-846b-d733daef1952",
   "metadata": {},
   "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",
   "id": "9817b380-8fd5-464e-84d9-29505ae7fd41",
   "metadata": {},
   "source": [
    "## Perform pycytominer pipeline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "cc03cbc9-3c0d-4251-a66c-a82a2f10c166",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Now performing pycytominer pipeline for Plate_3\n",
      "HET cells have been removed from Plate_3\n",
      "Now performing pycytominer pipeline for Plate_3_prime\n",
      "HET cells have been removed from Plate_3_prime\n"
     ]
    }
   ],
   "source": [
    "for plate, info in plate_info_dictionary.items():\n",
    "    print(f\"Now performing pycytominer pipeline for {plate}\")\n",
    "    output_aggregated_file = str(pathlib.Path(f\"{output_dir}/{plate}_bulk.parquet\"))\n",
    "    output_annotated_file = str(\n",
    "        pathlib.Path(f\"{output_dir}/{plate}_bulk_annotated.parquet\")\n",
    "    )\n",
    "    output_normalized_file = str(\n",
    "        pathlib.Path(f\"{output_dir}/{plate}_bulk_normalized.parquet\")\n",
    "    )\n",
    "    output_feature_select_file = str(\n",
    "        pathlib.Path(f\"{output_dir}/{plate}_bulk_feature_selected.parquet\")\n",
    "    )\n",
    "\n",
    "    # Load single-cell profiles\n",
    "    single_cell_df = pd.read_parquet(info[\"dest_path\"])\n",
    "\n",
    "    # Load platemap\n",
    "    platemap_df = pd.read_csv(info[\"platemap_path\"])\n",
    "\n",
    "    # Step 1: Aggregation\n",
    "    aggregate(\n",
    "        population_df=single_cell_df,\n",
    "        operation=\"median\",\n",
    "        strata=[\"Image_Metadata_Plate\", \"Image_Metadata_Well\"],\n",
    "        output_file=output_aggregated_file,\n",
    "        output_type=\"parquet\",\n",
    "    )\n",
    "\n",
    "    # Step 2: Annotation\n",
    "    annotated_df = annotate(\n",
    "        profiles=output_aggregated_file,\n",
    "        platemap=platemap_df,\n",
    "        join_on=[\"Metadata_well_position\", \"Image_Metadata_Well\"],\n",
    "    )\n",
    "\n",
    "    # For only plates 3 and 3 prime, remove any rows with HET due to contamination\n",
    "    if plate in [\"Plate_3\", \"Plate_3_prime\"]:\n",
    "        # Filter single-cell profiles, removing HET genotype\n",
    "        annotated_df = annotated_df[annotated_df[\"Metadata_genotype\"] != \"HET\"]\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",
    "    # Step 3: Normalization\n",
    "    normalized_df = normalize(\n",
    "        profiles=annotated_df,\n",
    "        method=\"standardize\",\n",
    "        output_file=output_normalized_file,\n",
    "        output_type=\"parquet\",\n",
    "    )\n",
    "\n",
    "    # Step 4: 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",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "4979678b-d0a0-40f8-b3e3-3ec5734ad8a7",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(48, 966)\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Metadata_WellRow</th>\n",
       "      <th>Metadata_WellCol</th>\n",
       "      <th>Metadata_gene_name</th>\n",
       "      <th>Metadata_genotype</th>\n",
       "      <th>Metadata_seed_density</th>\n",
       "      <th>Metadata_Plate</th>\n",
       "      <th>Metadata_Well</th>\n",
       "      <th>Cytoplasm_AreaShape_BoundingBoxMaximum_X</th>\n",
       "      <th>Cytoplasm_AreaShape_Compactness</th>\n",
       "      <th>Cytoplasm_AreaShape_Eccentricity</th>\n",
       "      <th>...</th>\n",
       "      <th>Nuclei_Texture_InfoMeas2_RFP_3_01_256</th>\n",
       "      <th>Nuclei_Texture_InfoMeas2_RFP_3_02_256</th>\n",
       "      <th>Nuclei_Texture_InfoMeas2_RFP_3_03_256</th>\n",
       "      <th>Nuclei_Texture_InverseDifferenceMoment_DAPI_3_00_256</th>\n",
       "      <th>Nuclei_Texture_InverseDifferenceMoment_GFP_3_00_256</th>\n",
       "      <th>Nuclei_Texture_InverseDifferenceMoment_RFP_3_00_256</th>\n",
       "      <th>Nuclei_Texture_InverseDifferenceMoment_RFP_3_01_256</th>\n",
       "      <th>Nuclei_Texture_InverseDifferenceMoment_RFP_3_02_256</th>\n",
       "      <th>Nuclei_Texture_InverseDifferenceMoment_RFP_3_03_256</th>\n",
       "      <th>Nuclei_Texture_SumEntropy_DAPI_3_00_256</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>B</td>\n",
       "      <td>1</td>\n",
       "      <td>NF1</td>\n",
       "      <td>WT</td>\n",
       "      <td>500</td>\n",
       "      <td>Plate_3_prime</td>\n",
       "      <td>B1</td>\n",
       "      <td>0.21263</td>\n",
       "      <td>-0.397138</td>\n",
       "      <td>0.223489</td>\n",
       "      <td>...</td>\n",
       "      <td>-1.545274</td>\n",
       "      <td>-1.031807</td>\n",
       "      <td>-0.253806</td>\n",
       "      <td>-0.391714</td>\n",
       "      <td>-0.345233</td>\n",
       "      <td>1.190598</td>\n",
       "      <td>0.629662</td>\n",
       "      <td>0.061780</td>\n",
       "      <td>2.371923</td>\n",
       "      <td>1.229115</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>C</td>\n",
       "      <td>1</td>\n",
       "      <td>NF1</td>\n",
       "      <td>WT</td>\n",
       "      <td>500</td>\n",
       "      <td>Plate_3_prime</td>\n",
       "      <td>C1</td>\n",
       "      <td>0.39603</td>\n",
       "      <td>0.658929</td>\n",
       "      <td>3.521152</td>\n",
       "      <td>...</td>\n",
       "      <td>-2.213429</td>\n",
       "      <td>-0.564723</td>\n",
       "      <td>-2.272657</td>\n",
       "      <td>-0.857327</td>\n",
       "      <td>-0.212326</td>\n",
       "      <td>-0.421005</td>\n",
       "      <td>0.259068</td>\n",
       "      <td>1.201338</td>\n",
       "      <td>1.318589</td>\n",
       "      <td>1.509185</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>2 rows × 966 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "  Metadata_WellRow  Metadata_WellCol Metadata_gene_name Metadata_genotype  \\\n",
       "0                B                 1                NF1                WT   \n",
       "1                C                 1                NF1                WT   \n",
       "\n",
       "   Metadata_seed_density Metadata_Plate Metadata_Well  \\\n",
       "0                    500  Plate_3_prime            B1   \n",
       "1                    500  Plate_3_prime            C1   \n",
       "\n",
       "   Cytoplasm_AreaShape_BoundingBoxMaximum_X  Cytoplasm_AreaShape_Compactness  \\\n",
       "0                                   0.21263                        -0.397138   \n",
       "1                                   0.39603                         0.658929   \n",
       "\n",
       "   Cytoplasm_AreaShape_Eccentricity  ...  \\\n",
       "0                          0.223489  ...   \n",
       "1                          3.521152  ...   \n",
       "\n",
       "   Nuclei_Texture_InfoMeas2_RFP_3_01_256  \\\n",
       "0                              -1.545274   \n",
       "1                              -2.213429   \n",
       "\n",
       "   Nuclei_Texture_InfoMeas2_RFP_3_02_256  \\\n",
       "0                              -1.031807   \n",
       "1                              -0.564723   \n",
       "\n",
       "   Nuclei_Texture_InfoMeas2_RFP_3_03_256  \\\n",
       "0                              -0.253806   \n",
       "1                              -2.272657   \n",
       "\n",
       "   Nuclei_Texture_InverseDifferenceMoment_DAPI_3_00_256  \\\n",
       "0                                          -0.391714      \n",
       "1                                          -0.857327      \n",
       "\n",
       "   Nuclei_Texture_InverseDifferenceMoment_GFP_3_00_256  \\\n",
       "0                                          -0.345233     \n",
       "1                                          -0.212326     \n",
       "\n",
       "   Nuclei_Texture_InverseDifferenceMoment_RFP_3_00_256  \\\n",
       "0                                           1.190598     \n",
       "1                                          -0.421005     \n",
       "\n",
       "   Nuclei_Texture_InverseDifferenceMoment_RFP_3_01_256  \\\n",
       "0                                           0.629662     \n",
       "1                                           0.259068     \n",
       "\n",
       "   Nuclei_Texture_InverseDifferenceMoment_RFP_3_02_256  \\\n",
       "0                                           0.061780     \n",
       "1                                           1.201338     \n",
       "\n",
       "   Nuclei_Texture_InverseDifferenceMoment_RFP_3_03_256  \\\n",
       "0                                           2.371923     \n",
       "1                                           1.318589     \n",
       "\n",
       "   Nuclei_Texture_SumEntropy_DAPI_3_00_256  \n",
       "0                                 1.229115  \n",
       "1                                 1.509185  \n",
       "\n",
       "[2 rows x 966 columns]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Check output file\n",
    "test_df = load_profiles(output_feature_select_file)\n",
    "\n",
    "print(test_df.shape)\n",
    "test_df.head(2)"
   ]
  }
 ],
 "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.8.19"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}