main/COMO.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "source": [
    "# COMO: Constraint-based Optomization of Metabolic Objectives\n",
    "\n",
    "COMO is used to build computational models that simulate the biochemical and phisiological processes that occur in a cell or organism, known as constraint-based metabolic models. The basic idea behind a constraint-based metabolic model is to use a set of constraints to place boundaries on the system being modeled. These constraints may include (but are not limited to) limits on the availability of nutrients, energy requirements, and the maximum rates of metabolic reactions. COMO imposes these constraints within a specific context. This context includes the cell or tissue type being modeled, along with its disease state. In addition to creating metabolic models, COMO serves as a platform to identify (1) drug targets and (2) repurposable drugs for metabolism-impacting diseases.\n",
    "\n",
    "\n",
    "This pipeline has everything necessary to build a model from any combination of the following sources:\n",
    "- Bulk RNA-seq (total and mRNA)\n",
    "- Single-cell RNA-seq\n",
    "- Proteomics\n",
    "\n",
    "\n",
    "COMO does not require programming experience to create models. However, every step of the pipeline is easily accessable to promote modification, addition, or replacement of analysis steps. In addition, this docker container comes pre-loaded with popular R and Python libraries; if you would like to use a library and cannot install it for any reason, please [request it on our GitHub page](https://github.com/HelikarLab/COMO)!\n",
    "\n",
    "\n",
    "<h2>\n",
    "<font color='red'>⚠️ WARNING ⚠️</font>\n",
    "</h2>\n",
    "\n",
    "If you terminate your session after running Docker, any changes you make *will <ins>**not**</ins> be saved*. Please mount a local directory to the docker image, [as instructed on the GitHub page](https://helikarlab.github.io/COMO/#choosing-a-tag), to prevent data loss.\n",
    "\n",
    "# Before Starting\n",
    "## Input Files\n",
    "The proper input files, dependent on the types of data you are using, must be loaded before model creation. Some example files are included to build metabolic models of naive, Th1, Th2, and Th17 T-cell subtypes, and identify targets for rheumatoid arthritis.\n",
    "\n",
    "### RNA-seq\n",
    "A correctly formatted folder named \"COMO inputs\" in the data directory. Proper inputs can be generated using our Snakemake pipeline, [FastqToGeneCounts](https://github.com/HelikarLab/FastqToGeneCounts), which is specifically designed for use with COMO. RNA sequencing data can be single-cell, or bulk, but the provided Snakemake pipeline does not process single-cell data as of now. If you are processing RNA-seq data with an alternate procedure or importing a pre-made gene count matrix, follow the instructions [listed under Step 1](#Importing-a-Pre-Generated-Counts-Matrix)\n",
    "\n",
    "### Proteomics\n",
    "A matrix of measurement values, where rows are protein names in Entrez format and columns are sample names\n",
    "\n",
    "## Configuration Information\n",
    "You should upload configuration files (in Excel format, `.xlsx`) to `data/config_sheets`. The sheet names in these configuration files should correspond to the context (tissue name, cell name, etc.). The data in each sheet contains the sample names to include in that context-specific model. These sample names should correspond to the column name in the source data matrix, which will be output (or uploaded, if you have your own data) to `data/data_matrices/MODEL-NAME`\n",
    "\n",
    "# Drug Target Identification\n",
    "\n",
    "1. Preprocess Bulk RNA-seq data\n",
    "    1. Convert STAR-output gene count files into a unified matrix\n",
    "    2. Fetch necessary information about each gene in the matrix\n",
    "    3. Generate a configuration file\n",
    "2. Analyze any combination of RNA-seq or proteomics data, and output a list of active genes for each strategy\n",
    "3. Check for a consensus amongst strategies according to a desired rigor and merge into a singular set of active genes\n",
    "4. Create a tissue-specific model based on the list of active genes (from Step 3)\n",
    "5. Identify differential gene expression from disease datasets using RNA-seq transcriptomics information\n",
    "6. Identify drug targets and repurposable drugs. This step consists of four substeps:\n",
    "    1. Map drugs to models\n",
    "    2. Knock-out simulation\n",
    "    3. Compare results between perturbed and unperturbed models (i.e., knocked-out models vs non-knocked-out models)\n",
    "    4. Integrate with disease genes and create a score of drug targets"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "source": [
    "# Step 1: Data Preprocessing and Analysis\n",
    "\n",
    "The first step of COMO will perform processing and analysis on each of the following data:\n",
    "- Total RNA sequencing\n",
    "- mRNA sequencing\n",
    "- Proteomics\n",
    "\n",
    "## RNA-seq Data\n",
    "RNA sequencing data is read by COMO as a count matrix, where each column is a different sample or replicate named \"tissueName_SXRYrZ\", where:\n",
    "- \"`X`\" represents the study (or batch) number. Each study represents a new experiment\n",
    "- \"`Y`\" represents the replicate number\n",
    "- \"`Z`\" represents the run number. If the replicate does not contain multiple runs for a single replicate, then \"`rZ`\" should not be included.\n",
    "- \"`tissueName`\" represents the name of the model that will be built from this data. It should be consistent with other data sources if you would like them to be integrated.\n",
    "\n",
    "❗The `tissueName` identifier should not contain any special characters, including `_`. Doing so may interfere with parsing throughout this pipeline.\n",
    "\n",
    "Replicates should come from the same study or batch group. Different studies/batches can come from different published studies, as long as the tissue/cell was under similar enough conditions for your personal modeling purposes. \"Run numbers\" in the same replicate will be summed together.\n",
    "\n",
    "### Example\n",
    "Pretend `S1` represents a study done by Margaret and `S2` represents a different study done by John. Margaret's experiment contains three replicates, while John's only contains two. Each of these studies comes from m0 Macrophages. Using this cell name, we will set our tissue name to `m0Macro`. The studies were conducted in different labs, by different researches, at different points in time, even using different preparation kits. . Using this information, we have the following samples:\n",
    "\n",
    "<table style=\"border: 1px solid black; border-collapse: collapse;\">\n",
    "    <thead>\n",
    "        <tr>\n",
    "            <th colspan=\"1000\" style=\"text-align: center;\">m0 Macrophage Data</th>\n",
    "        </tr>\n",
    "    </thead>\n",
    "    <tbody>\n",
    "        <tr>\n",
    "            <td colspan=\"3\" style=\"padding: 10px; text-align: center; border-bottom: 1px solid black;\">Margaret's Data</td>\n",
    "            <td colspan=\"3\" style=\"padding: 10px; text-align: center; border-left: 1px solid black; border-bottom: 1px solid black;\">John's Data</td>\n",
    "        </tr>\n",
    "        <tr>\n",
    "            <td style=\"padding: 10px; text-align: center;\">Study</td>\n",
    "            <td style=\"padding: 10px; text-align: center;\">Replicate</td>\n",
    "            <td style=\"padding: 10px; text-align: center;\">Resulting Name</td>\n",
    "            <td style=\"padding: 10px; text-align: center; border-left: 1px solid black;\">Study</td>\n",
    "            <td style=\"padding: 10px; text-align: center;\">Replicate</td>\n",
    "            <td style=\"padding: 10px; text-align: center;\">Resulting Name</td>\n",
    "        </tr>\n",
    "        <tr>\n",
    "            <td style=\"padding: 10px; text-align: center;\">S1</td>\n",
    "            <td style=\"padding: 10px; text-align: center;\">R1</td>\n",
    "            <td style=\"padding: 10px; text-align: center;\">m0Macro_S1R1</td>\n",
    "            <td style=\"padding: 10px; text-align: center; border-left: 1px solid black;\">S2</td>\n",
    "            <td style=\"padding: 10px; text-align: center;\">R1</td>\n",
    "            <td style=\"padding: 10px; text-align: center;\">m0Macro_S2R1</td>\n",
    "        </tr>\n",
    "        <tr>\n",
    "            <td style=\"padding: 10px; text-align: center;\">S1</td>\n",
    "            <td style=\"padding: 10px; text-align: center;\">R2</td>\n",
    "            <td style=\"padding: 10px; text-align: center;\">m0Macro_S1R2</td>\n",
    "            <td style=\"padding: 10px; text-align: center; border-left: 1px solid black;\">S2</td>\n",
    "            <td style=\"padding: 10px; text-align: center;\">R2</td>\n",
    "            <td style=\"padding: 10px; text-align: center;\">m0Macro_S2R2</td>\n",
    "        </tr>\n",
    "        <tr>\n",
    "            <td style=\"padding: 10px; text-align: center;\">S1</td>\n",
    "            <td style=\"padding: 10px; text-align: center;\">R3</td>\n",
    "            <td style=\"padding: 10px; text-align: center;\">m0Macro_S1R3</td>\n",
    "            <td style=\"padding: 10px; text-align: center; border-left: 1px solid black;\">-</td>\n",
    "            <td style=\"padding: 10px; text-align: center;\">-</td>\n",
    "            <td style=\"padding: 10px; text-align: center;\">-</td>\n",
    "        </tr>\n",
    "    </tbody>\n",
    "</table>\n",
    "\n",
    "From the `Resulting Name` column, the `m0Macro_S1R1`, `m0Macro_S1R2`, and `m0Macro_S1R3` samples (Margaret's data) will be checked for gene expression consensus to generate a list of active genes in all three replicates. The same will be done for `m0Macro_S2R1` and `m0Macro_S2R2` (John's data). Once these two *separate* lists of active genes have been generated, expression *between* lists will be checked for additional consensus between the studies. This system is used not only to help maintain organization throughout COMO, but because most types of normalized gene counts cannot undergo direct comparisons across replicates. This is especially true for comparisons between different experiments. Therefore, COMO will convert normalized gene counts into a boolean list of active genes. These lists will be compared at the level of replicates in a study, and then again at the level of all provided studies. Finally, the active genes will be merged with the outputs of proteomics and various RNA-sequencing strategies if provided. The rigor used at each level is easily modifiable.\n",
    "\n",
    "\n",
    "### Initializing RNA-seq Data\n",
    "\n",
    "Please choose an option below:\n",
    "1. Importing a `COMO inputs` directory\n",
    "    1. [Initialization using the Snakemake Pipeline](https://github.com/HelikarLab/FastqToGeneCounts)\n",
    "    2. [Creating your own Inputs](#Creating-a-Properly-Formatted-COMO-inputs-Folder)\n",
    "2. [Importing a pre-generated gene counts file](#Importing-a-Pre-Generated-Counts-Matrix)\n",
    "\n",
    "#### Snakemake Pipeline\n",
    "It is recommended you use the available Snakemake pipeline to align to create a properly formatted `COMO inputs` folder. The pipeline also runs a series of quality control steps to help determine if any of the provided samples are not suitable for model creation. This pipeline can be found at https://github.com/HelikarLab/FastqToGeneCounts.\n",
    "\n",
    "The folder output from the snakemake pipeline can be uploaded directly to the folder `data/COMO inputs` in this pipeline\n",
    "\n",
    "Once this is done, continue to the code block at the end of this section\n",
    "\n",
    "#### Creating a Properly Formatted `COMO inputs` Folder\n",
    "\n",
    "\n",
    "If you are using your own alignment protocol, follow this section to create a properly formatted `COMO inputs` folder.\n",
    "\n",
    "The top-level of the directory will have separate tissue/cell types that models should be created from. The next level must have a folder called `geneCounts`, and optionally a `strandedness` folder. If you are using zFPKM normalization, two additional folders must be included: `layouts` and `fragmentSizes`. Inside each of these folders should be folders named `SX`, where `X` is a number that replicates are associated with.\n",
    "\n",
    "<br>\n",
    "\n",
    "<ins>Gene Counts</ins>\n",
    "Create a folder named `geneCounts`. The outputs of the STAR aligner using the `-quantMode GeneCounts` option should be included inside the \"study-number\" folders (`SX`) of `geneCounts`. To help you (and COMO!) stay organized, these outputs should be renamed `tissueName_SXRYrZ.tab`. Just like above, `X` is the study number, `Y` is the replicate number, and (if present), `Z` is the run number. If the replicate does not contain multiple runs, the `rZ` should be excluded from the name. Replicates should come from the same study/sample group. Different samples can come from different published studies as long as the experiments were performed under similar enough conditions for your modeling purposes.\n",
    "\n",
    "<ins>Strandedness</ins>\n",
    "Create a folder named `strandedness`. This folder should contain files named `tissueName_SXRYrZ_strandedness.txt`. These files must tell the strandedness of the RNA-sequencing method used. It should contain one of the following texts (and nothing else):\n",
    "    - `NONE`: If you don't know the strandedness\n",
    "    - `FIRST_READ_TRANSCRIPTION_STRAND`: If this RNA-sequencing sample originates from the first strand of cDNA, or the \"antisense\" strand\n",
    "    - `SECOND_READ_TRANSCRIPTION_STRAND`: If this RNA-sequencing sample originates from the second strand of cDNA, or the \"sense\" strand\n",
    "\n",
    "<ins>Layouts</ins>\n",
    "Create a folder a folder named `layouts`. Files should be named `tissueName_SXRYrZ_layout.txt, where each file tells the layout of the library used. It must contain one of the following texts, and nothing else:\n",
    "- `paired-end`: Paired-end reads were generated\n",
    "- `single-end`: Single-end reads were generated\n",
    "\n",
    "<ins>Fragment Sizes</ins>\n",
    "Create a folder named `fragmentSizes`. Files should be named `tissueName_SXRYrZ_fragment_sizes.txt` and contain the output of [RSeQC](https://rseqc.sourceforge.net/)'s `como/RNA_fragment_size.py` function.\n",
    "\n",
    "<ins>Preparation Methods</ins>\n",
    "Create a folder named `prepMethods`. Files should be named `tissueName_SXRYrZ_prep_method.txt`. Each file should tell the library preparation strategy. It must contain one of the following texts, and nothing else:\n",
    "- `total`: All mRNA expression was measured (mRNA, ncRNA, rRNA, etc.)\n",
    "- `mRNA`: Only polyA mRNA expression was measured\n",
    "\n",
    "It should be noted that these strategies only serve to differentiate the methods in the event that both are used to build a model. If a different library strategy is desired, you have two options:\n",
    "1. Replace one of these with a placeholder. If you only have polyA mRNA expression, you only have to enter data for those samples. Do not fill out any samples with `total`.\n",
    "2. With a little Python knowledge, a new strategy can easily be added to the `como/merge_xomics.py` file. If you would like to do so, the file is located under `como/merge_xomics.py` in this Jupyter Notebook\n",
    "\n",
    "#### Importing a Pre-Generated Counts Matrix\n",
    "Import a properly formatted counts matrix to `data/data_matrices/exampleTissue/gene_counts_matrix_exampleTissue.csv`. The rows should be named `exampleTissue_SXRY` (note the lack of a run number (`rZ`), runs should be summed into each replicate). If you are providing the count matrix this way, instead of generating one using the snakemake pipeline mentioned above, you must create a configuration file that has each sample's name, study number, and if using zFPKM, layout and mean fragment length. Use the provided template below to create yours. Once you have created this file and placed it under the `data/data_matrices/exampleTissue` directory, run the `como/rnaseq_preprocess.py` file with `preprocess-mode` set to `provide-matrix`.\n",
    "\n",
    "This method is best if you are downloading a premade count matrix, or using single-cell data that has already been batch corrected, clustered, and sorted into only the cell type of interest!\n",
    "\n",
    "\n",
    "<table style=\"border: 1px solid black; border-collapse: collapse;\">\n",
    "    <thead>\n",
    "        <tr>\n",
    "            <th colspan=\"1000\" style=\"text-align: center; border-bottom: 1px solid black;\">Example Gene Count Table</th>\n",
    "        </tr>\n",
    "    </thead>\n",
    "    <tbody>\n",
    "        <tr>\n",
    "            <td style=\"padding: 10px; text-align: center; border-bottom: 1px solid black;\">genes</td>\n",
    "            <td style=\"padding: 10px; text-align: center; border-left: 1px solid black; border-bottom: 1px solid black;\">exampleTissue_S1R1</td>\n",
    "            <td style=\"padding: 10px; text-align: center; border-left: 1px solid black; border-bottom: 1px solid black;\">exampleTissue_S1R2</td>\n",
    "            <td style=\"padding: 10px; text-align: center; border-left: 1px solid black; border-bottom: 1px solid black;\">exampleTissue_S2R1</td>\n",
    "            <td style=\"padding: 10px; text-align: center; border-left: 1px solid black; border-bottom: 1px solid black;\">exampleTissue_S2R2</td>\n",
    "        </tr>\n",
    "        <tr>\n",
    "            <td style=\"padding: 10px; text-align: center;\">ENSG00000000003</td>\n",
    "            <td style=\"padding: 10px; text-align: center;\">20</td>\n",
    "            <td style=\"padding: 10px; text-align: center;\">29</td>\n",
    "            <td style=\"padding: 10px; text-align: center;\">52</td>\n",
    "            <td style=\"padding: 10px; text-align: center;\">71</td>\n",
    "        </tr>\n",
    "        <tr>\n",
    "            <td style=\"padding: 10px; text-align: center;\">ENSG00000000005</td>\n",
    "            <td style=\"padding: 10px; text-align: center;\">0</td>\n",
    "            <td style=\"padding: 10px; text-align: center;\">0</td>\n",
    "            <td style=\"padding: 10px; text-align: center;\">0</td>\n",
    "            <td style=\"padding: 10px; text-align: center;\">0</td>\n",
    "        </tr>\n",
    "        <tr>\n",
    "            <td style=\"padding: 10px; text-align: center;\">ENSG00000000419</td>\n",
    "            <td style=\"padding: 10px; text-align: center;\">1354</td>\n",
    "            <td style=\"padding: 10px; text-align: center;\">2081</td>\n",
    "            <td style=\"padding: 10px; text-align: center;\">1760</td>\n",
    "            <td style=\"padding: 10px; text-align: center;\">3400</td>\n",
    "        </tr>\n",
    "    </tbody>\n",
    "</table>"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-12-16T20:09:29.761117Z",
     "start_time": "2024-12-16T20:09:29.754605Z"
    }
   },
   "source": [
    "from pathlib import Path\n",
    "\n",
    "from como.data_types import RNAType\n",
    "from como.rnaseq_preprocess import rnaseq_preprocess\n",
    "\n",
    "taxon_id = 9606\n",
    "context_names = [\"naiveB\"]\n",
    "\n",
    "# Input Components\n",
    "como_context_dir = {context: Path(f\"data/COMO_input/{context}\") for context in context_names}\n",
    "gene_info_filepath = {context: Path(f\"data/results/{context}/gene_info.csv\") for context in context_names}\n",
    "trna_matrix_filepath = {context: Path(f\"data/results/{context}/total-rna/totalrna_{context}.csv\") for context in context_names}\n",
    "mrna_matrix_filepath = {context: Path(f\"data/results/{context}/mrna/mrna_{context}.csv\") for context in context_names}\n",
    "proteomics_matrix_filepath = {context: Path(f\"data/data_matrices/{context}/protein_abundance_{context}.csv\") for context in context_names}\n",
    "\n",
    "# No single-cell data is provided by default; COMO accepts single-cell data in CSV or h5ad format\n",
    "# If you are using single-cell data, adjust the following lines accordingly\n",
    "scrna_matrix_filepath = {context: Path(f\"data/results/{context}/scrna/scrna_{context}.csv\") for context in context_names}\n",
    "# scrna_matrix_filepath = [Path(f\"data/results/{context}/scrna/scrna_{context}.h5ad\") for context in context_names]\n",
    "\n",
    "trna_metadata_filepath = Path(\"data/config_sheets/trna_config.xlsx\")\n",
    "mrna_metadata_filepath = Path(\"data/config_sheets/mrna_config.xlsx\")\n",
    "proteomics_metadata_filepath = Path(\"data/config_sheets/proteomics_config.xlsx\")\n",
    "\n",
    "# Output Components\n"
   ],
   "outputs": [],
   "execution_count": 11
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### RNA-seq Preprocessing Parameters\n",
    "- `context_names`: The tissue/cell types to use. This is a simple space-separated list of items, such as \"naiveB regulatoryTcell\"\n",
    "- `gene_format`: The format of input genes, accepts `\"Extrez\"`, `\"Emsembl\"` or `\"Symbol\"`\n",
    "- `taxon_id`: The [NCBI Taxon ID](https://www.ncbi.nlm.nih.gov/taxonomy) to use\n",
    "- `preprocess_mode`: This should be set to `\"create-matrix\"` if you are **not** providing a matrix, otherwise set it to `\"provide-matrix\"`"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-12-16T20:07:08.331063Z",
     "start_time": "2024-12-16T20:06:10.188429Z"
    }
   },
   "source": [
    "for context in context_names:\n",
    "    if context not in {*trna_matrix_filepath, *mrna_matrix_filepath}:\n",
    "        continue\n",
    "    await rnaseq_preprocess(\n",
    "        context_name=context,\n",
    "        taxon=taxon_id,\n",
    "        como_context_dir=como_context_dir[context],\n",
    "        input_matrix_filepath=None,\n",
    "        output_gene_info_filepath=gene_info_filepath[context],\n",
    "        output_trna_metadata_filepath=trna_metadata_filepath,\n",
    "        output_trna_count_matrix_filepath=trna_matrix_filepath[context],\n",
    "        output_mrna_metadata_filepath=mrna_metadata_filepath,\n",
    "        output_mrna_count_matrix_filepath=mrna_matrix_filepath[context],\n",
    "        cache=True,\n",
    "        log_level=\"INFO\",\n",
    "    )"
   ],
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\u001B[32m2024-12-16 14:06:13\u001B[0m | \u001B[32m\u001B[1mSUCCESS \u001B[0m | \u001B[36mcomo.rnaseq_preprocess\u001B[0m:\u001B[36m277\u001B[0m - \u001B[32m\u001B[1mWrote gene count matrix for 'mrna' RNA at '/Users/joshl/Projects/COMO/main/data/results/naiveB/polya-rna/polyarna_naiveB.csv'\u001B[0m\n",
      "\u001B[32m2024-12-16 14:06:13\u001B[0m | \u001B[32m\u001B[1mSUCCESS \u001B[0m | \u001B[36mcomo.rnaseq_preprocess\u001B[0m:\u001B[36m277\u001B[0m - \u001B[32m\u001B[1mWrote gene count matrix for 'total' RNA at '/Users/joshl/Projects/COMO/main/data/results/naiveB/total-rna/totalrna_naiveB.csv'\u001B[0m\n",
      "\u001B[32m2024-12-16 14:06:13\u001B[0m | \u001B[1mINFO    \u001B[0m | \u001B[36mcomo.rnaseq_preprocess\u001B[0m:\u001B[36m468\u001B[0m - \u001B[1mFetching gene info (this may take 1-5 minutes depending on the number of genes and your internet connection)\u001B[0m\n",
      "\u001B[32m2024-12-16 14:07:08\u001B[0m | \u001B[32m\u001B[1mSUCCESS \u001B[0m | \u001B[36mcomo.rnaseq_preprocess\u001B[0m:\u001B[36m505\u001B[0m - \u001B[32m\u001B[1mGene Info file written at '/Users/joshl/Projects/COMO/main/data/results/naiveB/gene_info.csv'\u001B[0m\n"
     ]
    }
   ],
   "execution_count": 8
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## Identification of Gene Activity in Transcriptomic and Proteomic Datasets\n",
    "\n",
    "This part of Step 1 will identify gene activity in the following data sources:\n",
    "- RNA-seq (total, mRNA, or single cell)\n",
    "- Proteomics\n",
    "\n",
    "Only one source is required for model generation, but multiple sources can be helpful for additional validation if they are of high enough quality\n",
    "\n",
    "### Filtering Raw Counts\n",
    "Regardless of normalization technique used, or provided files used for RNA-seq, preprocessing is required to fetch relevent gene information needed for harmonization and normalization such as Entrez ID, and the start and end postions. Currently, COMO can filter raw RNA-sequencing counts using one of the following normalization techniques:\n",
    "\n",
    "#### Transcripts Per Million Quantile\n",
    "TPM Quantile. Each replicate is normalized with Transcripts-Per-Million, and an upper quantiile is taken to create a boolean list of active genes for the replicate (i.e., `R1`). Replicates are compared for consensus within the study, and then studies are compared between one another for additional consensus. The strictness of the consensus easily be set using the appropriate option within the `rnaseq_gen.py` code-block.\n",
    "\n",
    "This method is recommended if you want more control over the size of the model; smaller models can include only the most expressed reactions, and larger models can encompass less essentail reactions\n",
    "\n",
    "#### zFPKM\n",
    "This method is outlined by [Hart et. al](https://pubmed.ncbi.nlm.nih.gov/24215113/). Counts will be normalized using zFPKM and genes > -3 will be considered \"expressed\" per Hart's recommendation. Expressed genes will be checked for consensus at the replicate and study level.\n",
    "\n",
    "This method is recommended if you want to less control over which genes are essential, and instead use the most standardized method of active gene determination. This method is more \"hands-off\" than the above TPM Quantile method.\n",
    "\n",
    "#### Counts Per Million\n",
    "This is a flat cutoff value of counts per million normalized values. Gene expression will be checked for consensus at the replicate and study level.\n",
    "\n",
    "This method is not recommended, as zFPKM is much more robust for a similar level of \"hands-off\" model building\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### RNA Sequencing Analysis\n",
    "#### Bulk RNA Sequencing\n",
    "\n",
    "This has multiple strategies of library preparation (total, polyA-mRNA). If you are using public data, you may encounter a situation where you would like to use a combination of bulk RNA sequencing data produced using two different library preparation strategies.\n",
    "\n",
    "COMO currently supports the two most common strategies, mRNA polyA enriched RNA sequencing, and total RNA sequencing. Because of the expected differences in distribution of transcripts, COMO is written to handle each strategy seperately before the integration step. The recommended Snakemake alignment pipeline is designed to work with COMO's preprocessing step ([Step 1, above](Step-1:-Initialize-and-Preprocess-RNA-seq-data)) to split RNA sequencing data from GEO into seperate input matrices and configuration files.\n",
    "\n",
    "To create a gene expression file for total RNA sequencing data, use `\"total\"` for the \"`--library-prep`\" argument.\n",
    "To create a gene expression file for mRNA polyA enriched data, use `\"mRNA\"` for the  \"`--library-prep`\" argument.\n",
    "\n",
    "The analysis of each strategy is identical. Specifying the type of analysis (total vs mRNA) only serves to ensure COMO analyzes them seperately.\n",
    "\n",
    "#### Single Cell RNA Sequencing\n",
    "While the Snakemake pipeline does not yet support single-cell alignment, and COMO does not yet support automated configuration file and counts matrix file creation for single-cell alignment output from STAR, it is possible to use single-cell data to create a model with COMO. Because normalization strategies can be applied to single-cell data in the same way it is applied to bulk RNA sequencing, `como/rnaseq_gen.py` can be used with a provided counts matrix and configuration file, from [Step 1](Step-1:-Initialize-and-Preprocess-RNA-seq-data), above. Just like `\"total\"` and `\"mRNA\"`, `como/rnaseq_gen.py` can be executed with `\"SC\"` as the \"`--library-prep`\" argument to help COMO differentiate it from any bulk RNA sequencing data if multiple strategies are being used."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "### Total RNA Sequencing Generation\n",
    "#### Parameters\n",
    "- `trnaseq_config_file`: The configuration filename for total RNA. This file is found under the `data/config_sheets` folder\n",
    "- `rep_ratio`: The proportion of replicates before a gene is considered \"active\" in a study\n",
    "- `group_ratio`: The proportion of studies with expression required for a gene to be considered \"active\"\n",
    "- `rep_ratio_h`: The proportion of replicates that must express a gene before that gene is considered \"high-confidience\"\n",
    "- `group_ratio_h`: The proportion of studies that must express a gene before that gene is considered \"high-confidence\"\n",
    "- `technique`: The technique to use. Options are: `\"quantile\"`, `\"cpm\"`, or `\"zfpkm\"`. The difference in these options is discussed above\n",
    "- `quantile`: The cutoff Transcripts-Per-Million quantile for filtering\n",
    "- `min_zfpkm`: The cutoff for Counts-Per-Million filtering\n",
    "- `prep_method`: The library method used for preparation. Options are: `\"total\"`, `\"mRNA\"`, or `\"SC\"`,\n"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-12-16T20:08:00.274368Z",
     "start_time": "2024-12-16T20:07:48.637980Z"
    }
   },
   "source": [
    "from como.rnaseq_gen import FilteringTechnique, rnaseq_gen\n",
    "\n",
    "replicate_ratio = 0.75\n",
    "high_confidence_replicate_ratio = 1.0\n",
    "batch_ratio = 0.75\n",
    "high_confidence_batch_ratio = 1.0\n",
    "technique = FilteringTechnique.ZFPKM\n",
    "cutoff = -3\n",
    "\n",
    "\n",
    "for context in context_names:\n",
    "    if context not in trna_matrix_filepath:\n",
    "        continue\n",
    "    await rnaseq_gen(\n",
    "        context_name=context,\n",
    "        input_rnaseq_filepath=trna_matrix_filepath[context],\n",
    "        input_gene_info_filepath=gene_info_filepath[context],\n",
    "        output_boolean_activity_filepath=trna_matrix_filepath[context],\n",
    "        prep=RNAType.TRNA,\n",
    "        taxon_id=taxon_id,\n",
    "        input_metadata_filepath=trna_metadata_filepath,\n",
    "        replicate_ratio=replicate_ratio,\n",
    "        high_replicate_ratio=high_confidence_replicate_ratio,\n",
    "        batch_ratio=batch_ratio,\n",
    "        high_batch_ratio=high_confidence_batch_ratio,\n",
    "        technique=technique,\n",
    "        cutoff=cutoff\n",
    "    )"
   ],
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/joshl/Projects/COMO/main/como/rnaseq_gen.py:430: FutureWarning: The behavior of DataFrame concatenation with empty or all-NA entries is deprecated. In a future version, this will no longer exclude empty or all-NA columns when determining the result dtypes. To retain the old behavior, exclude the relevant entries before the concat operation.\n",
      "  mega_df = pd.concat([mega_df, df], ignore_index=True)\n",
      "\u001B[32m2024-12-16 14:08:00\u001B[0m | \u001B[1mINFO    \u001B[0m | \u001B[36mcomo.rnaseq_gen\u001B[0m:\u001B[36m691\u001B[0m - \u001B[1mnaiveB - Found 9492 expressed and 9492 confidently expressed genes\u001B[0m\n",
      "\u001B[32m2024-12-16 14:08:00\u001B[0m | \u001B[32m\u001B[1mSUCCESS \u001B[0m | \u001B[36mcomo.rnaseq_gen\u001B[0m:\u001B[36m694\u001B[0m - \u001B[32m\u001B[1mWrote boolean matrix to data/results/naiveB/total-rna/totalrna_naiveB.csv\u001B[0m\n"
     ]
    }
   ],
   "execution_count": 9
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## mRNA Sequencing Generation\n",
    "These parameters are identical to the ones listed for [total RNA sequencing](#Total-RNA-Sequencing-Generation), but they are listed again here for ease of reference\n",
    "\n",
    "### Parameters\n",
    "- `mrnaseq_config_file`: The configuration filename for total RNA. This file is found under the `data/config_sheets` folder\n",
    "- `rep_ratio`: The proportion of replicates before a gene is considered \"active\" in a study\n",
    "- `group_ratio`: The proportion of studies with expression required for a gene to be considered \"active\"\n",
    "- `rep_ratio_h`: The proportion of replicates that must express a gene before that gene is considered \"high-confidience\"\n",
    "- `group_ratio_h`: The proportion of studies that must express a gene before that gene is considered \"high-confidence\"\n",
    "- `technique`: The technique to use. Options are: `\"quantile\"`, `\"cpm\"`, or `\"zfpkm\"`. The difference in these options is discussed above\n",
    "- `quantile`: The cutoff Transcripts-Per-Million quantile for filtering\n",
    "- `min_zfpkm`: The cutoff for Counts-Per-Million filtering\n",
    "- `prep_method`: The library method used for preparation. Options are: `\"total\"`, `\"mRNA\"`, or `\"SC\"`,\n"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-12-16T20:15:16.139743Z",
     "start_time": "2024-12-16T20:15:15.998158Z"
    }
   },
   "source": [
    "from como.rnaseq_gen import FilteringTechnique, rnaseq_gen\n",
    "\n",
    "replicate_ratio = 0.75\n",
    "high_confidence_replicate_ratio = 1.0\n",
    "batch_ratio = 0.75\n",
    "high_confidence_batch_ratio = 1.0\n",
    "technique = FilteringTechnique.ZFPKM\n",
    "cutoff = -3\n",
    "\n",
    "for context in context_names:\n",
    "    if context not in mrna_matrix_filepath:\n",
    "        continue\n",
    "    await rnaseq_gen(\n",
    "        context_name=context,\n",
    "        input_rnaseq_filepath=mrna_matrix_filepath[context],\n",
    "        input_gene_info_filepath=gene_info_filepath[context],\n",
    "        output_boolean_activity_filepath=mrna_matrix_filepath[context],\n",
    "        prep=RNAType.MRNA,\n",
    "        taxon_id=taxon_id,\n",
    "        input_metadata_filepath=mrna_metadata_filepath,\n",
    "        replicate_ratio=replicate_ratio,\n",
    "        high_replicate_ratio=high_confidence_replicate_ratio,\n",
    "        batch_ratio=batch_ratio,\n",
    "        high_batch_ratio=high_confidence_batch_ratio,\n",
    "        technique=technique,\n",
    "        cutoff=cutoff\n",
    "    )"
   ],
   "outputs": [
    {
     "ename": "KeyError",
     "evalue": "'ensembl_gene_id'",
     "output_type": "error",
     "traceback": [
      "\u001B[0;31m---------------------------------------------------------------------------\u001B[0m",
      "\u001B[0;31mKeyError\u001B[0m                                  Traceback (most recent call last)",
      "File \u001B[0;32m~/Projects/COMO/.venv/lib/python3.10/site-packages/pandas/core/indexes/base.py:3805\u001B[0m, in \u001B[0;36mIndex.get_loc\u001B[0;34m(self, key)\u001B[0m\n\u001B[1;32m   3804\u001B[0m \u001B[38;5;28;01mtry\u001B[39;00m:\n\u001B[0;32m-> 3805\u001B[0m     \u001B[38;5;28;01mreturn\u001B[39;00m \u001B[38;5;28;43mself\u001B[39;49m\u001B[38;5;241;43m.\u001B[39;49m\u001B[43m_engine\u001B[49m\u001B[38;5;241;43m.\u001B[39;49m\u001B[43mget_loc\u001B[49m\u001B[43m(\u001B[49m\u001B[43mcasted_key\u001B[49m\u001B[43m)\u001B[49m\n\u001B[1;32m   3806\u001B[0m \u001B[38;5;28;01mexcept\u001B[39;00m \u001B[38;5;167;01mKeyError\u001B[39;00m \u001B[38;5;28;01mas\u001B[39;00m err:\n",
      "File \u001B[0;32mindex.pyx:167\u001B[0m, in \u001B[0;36mpandas._libs.index.IndexEngine.get_loc\u001B[0;34m()\u001B[0m\n",
      "File \u001B[0;32mindex.pyx:196\u001B[0m, in \u001B[0;36mpandas._libs.index.IndexEngine.get_loc\u001B[0;34m()\u001B[0m\n",
      "File \u001B[0;32mpandas/_libs/hashtable_class_helper.pxi:7081\u001B[0m, in \u001B[0;36mpandas._libs.hashtable.PyObjectHashTable.get_item\u001B[0;34m()\u001B[0m\n",
      "File \u001B[0;32mpandas/_libs/hashtable_class_helper.pxi:7089\u001B[0m, in \u001B[0;36mpandas._libs.hashtable.PyObjectHashTable.get_item\u001B[0;34m()\u001B[0m\n",
      "\u001B[0;31mKeyError\u001B[0m: 'ensembl_gene_id'",
      "\nThe above exception was the direct cause of the following exception:\n",
      "\u001B[0;31mKeyError\u001B[0m                                  Traceback (most recent call last)",
      "Cell \u001B[0;32mIn[15], line 13\u001B[0m\n\u001B[1;32m     11\u001B[0m \u001B[38;5;28;01mif\u001B[39;00m context \u001B[38;5;129;01mnot\u001B[39;00m \u001B[38;5;129;01min\u001B[39;00m mrna_matrix_filepath:\n\u001B[1;32m     12\u001B[0m     \u001B[38;5;28;01mcontinue\u001B[39;00m\n\u001B[0;32m---> 13\u001B[0m \u001B[38;5;28;01mawait\u001B[39;00m rnaseq_gen(\n\u001B[1;32m     14\u001B[0m     context_name\u001B[38;5;241m=\u001B[39mcontext,\n\u001B[1;32m     15\u001B[0m     input_rnaseq_filepath\u001B[38;5;241m=\u001B[39mmrna_matrix_filepath[context],\n\u001B[1;32m     16\u001B[0m     input_gene_info_filepath\u001B[38;5;241m=\u001B[39mgene_info_filepath[context],\n\u001B[1;32m     17\u001B[0m     output_rnaseq_filepath\u001B[38;5;241m=\u001B[39mmrna_matrix_filepath[context],\n\u001B[1;32m     18\u001B[0m     prep\u001B[38;5;241m=\u001B[39mRNAPrepMethod\u001B[38;5;241m.\u001B[39mMRNA,\n\u001B[1;32m     19\u001B[0m     taxon\u001B[38;5;241m=\u001B[39mtaxon_id,\n\u001B[1;32m     20\u001B[0m     input_metadata_filepath\u001B[38;5;241m=\u001B[39mmrna_metadata_filepath,\n\u001B[1;32m     21\u001B[0m     replicate_ratio\u001B[38;5;241m=\u001B[39mreplicate_ratio,\n\u001B[1;32m     22\u001B[0m     high_replicate_ratio\u001B[38;5;241m=\u001B[39mhigh_confidence_replicate_ratio,\n\u001B[1;32m     23\u001B[0m     batch_ratio\u001B[38;5;241m=\u001B[39mbatch_ratio,\n\u001B[1;32m     24\u001B[0m     high_batch_ratio\u001B[38;5;241m=\u001B[39mhigh_confidence_batch_ratio,\n\u001B[1;32m     25\u001B[0m     technique\u001B[38;5;241m=\u001B[39mtechnique,\n\u001B[1;32m     26\u001B[0m     cutoff\u001B[38;5;241m=\u001B[39mcutoff\n\u001B[1;32m     27\u001B[0m )\n",
      "File \u001B[0;32m~/Projects/COMO/main/como/rnaseq_gen.py:787\u001B[0m, in \u001B[0;36mrnaseq_gen\u001B[0;34m(context_name, input_rnaseq_filepath, input_gene_info_filepath, output_rnaseq_filepath, prep, taxon, input_metadata_filepath, input_metadata_df, replicate_ratio, high_replicate_ratio, batch_ratio, high_batch_ratio, technique, cutoff)\u001B[0m\n\u001B[1;32m    784\u001B[0m logger\u001B[38;5;241m.\u001B[39mdebug(\u001B[38;5;124mf\u001B[39m\u001B[38;5;124m\"\u001B[39m\u001B[38;5;124mStarting \u001B[39m\u001B[38;5;124m'\u001B[39m\u001B[38;5;132;01m{\u001B[39;00mcontext_name\u001B[38;5;132;01m}\u001B[39;00m\u001B[38;5;124m'\u001B[39m\u001B[38;5;124m\"\u001B[39m)\n\u001B[1;32m    785\u001B[0m output_rnaseq_filepath\u001B[38;5;241m.\u001B[39mparent\u001B[38;5;241m.\u001B[39mmkdir(parents\u001B[38;5;241m=\u001B[39m\u001B[38;5;28;01mTrue\u001B[39;00m, exist_ok\u001B[38;5;241m=\u001B[39m\u001B[38;5;28;01mTrue\u001B[39;00m)\n\u001B[0;32m--> 787\u001B[0m \u001B[38;5;28;01mawait\u001B[39;00m _save_rnaseq_tests(\n\u001B[1;32m    788\u001B[0m     context_name\u001B[38;5;241m=\u001B[39mcontext_name,\n\u001B[1;32m    789\u001B[0m     rnaseq_matrix\u001B[38;5;241m=\u001B[39m\u001B[38;5;28;01mawait\u001B[39;00m _read_counts(input_rnaseq_filepath),\n\u001B[1;32m    790\u001B[0m     metadata_df\u001B[38;5;241m=\u001B[39minput_metadata_df \u001B[38;5;129;01mor\u001B[39;00m \u001B[38;5;28;01mawait\u001B[39;00m _create_metadata_df(input_metadata_filepath),\n\u001B[1;32m    791\u001B[0m     gene_info_df\u001B[38;5;241m=\u001B[39mpd\u001B[38;5;241m.\u001B[39mread_csv(input_gene_info_filepath),\n\u001B[1;32m    792\u001B[0m     output_filepath\u001B[38;5;241m=\u001B[39moutput_rnaseq_filepath,\n\u001B[1;32m    793\u001B[0m     prep\u001B[38;5;241m=\u001B[39mprep,\n\u001B[1;32m    794\u001B[0m     taxon\u001B[38;5;241m=\u001B[39mtaxon,\n\u001B[1;32m    795\u001B[0m     replicate_ratio\u001B[38;5;241m=\u001B[39mreplicate_ratio,\n\u001B[1;32m    796\u001B[0m     batch_ratio\u001B[38;5;241m=\u001B[39mbatch_ratio,\n\u001B[1;32m    797\u001B[0m     high_replicate_ratio\u001B[38;5;241m=\u001B[39mhigh_replicate_ratio,\n\u001B[1;32m    798\u001B[0m     high_batch_ratio\u001B[38;5;241m=\u001B[39mhigh_batch_ratio,\n\u001B[1;32m    799\u001B[0m     technique\u001B[38;5;241m=\u001B[39mtechnique,\n\u001B[1;32m    800\u001B[0m     cut_off\u001B[38;5;241m=\u001B[39mcutoff,\n\u001B[1;32m    801\u001B[0m )\n",
      "File \u001B[0;32m~/Projects/COMO/main/como/rnaseq_gen.py:645\u001B[0m, in \u001B[0;36m_save_rnaseq_tests\u001B[0;34m(context_name, rnaseq_matrix, metadata_df, gene_info_df, output_filepath, prep, taxon, replicate_ratio, batch_ratio, high_replicate_ratio, high_batch_ratio, technique, cut_off)\u001B[0m\n\u001B[1;32m    636\u001B[0m \u001B[38;5;250m\u001B[39m\u001B[38;5;124;03m\"\"\"Save the results of the RNA-Seq tests to a CSV file.\"\"\"\u001B[39;00m\n\u001B[1;32m    637\u001B[0m filtering_options \u001B[38;5;241m=\u001B[39m _FilteringOptions(\n\u001B[1;32m    638\u001B[0m     replicate_ratio\u001B[38;5;241m=\u001B[39mreplicate_ratio,\n\u001B[1;32m    639\u001B[0m     batch_ratio\u001B[38;5;241m=\u001B[39mbatch_ratio,\n\u001B[0;32m   (...)\u001B[0m\n\u001B[1;32m    642\u001B[0m     high_batch_ratio\u001B[38;5;241m=\u001B[39mhigh_batch_ratio,\n\u001B[1;32m    643\u001B[0m )\n\u001B[0;32m--> 645\u001B[0m read_counts_results: _ReadMatrixResults \u001B[38;5;241m=\u001B[39m \u001B[38;5;28;01mawait\u001B[39;00m _build_matrix_results(\n\u001B[1;32m    646\u001B[0m     matrix\u001B[38;5;241m=\u001B[39mrnaseq_matrix,\n\u001B[1;32m    647\u001B[0m     gene_info\u001B[38;5;241m=\u001B[39mgene_info_df,\n\u001B[1;32m    648\u001B[0m     metadata_df\u001B[38;5;241m=\u001B[39mmetadata_df,\n\u001B[1;32m    649\u001B[0m     taxon\u001B[38;5;241m=\u001B[39mtaxon,\n\u001B[1;32m    650\u001B[0m )\n\u001B[1;32m    651\u001B[0m metrics \u001B[38;5;241m=\u001B[39m read_counts_results\u001B[38;5;241m.\u001B[39mmetrics\n\u001B[1;32m    652\u001B[0m entrez_gene_ids \u001B[38;5;241m=\u001B[39m read_counts_results\u001B[38;5;241m.\u001B[39mentrez_gene_ids\n",
      "File \u001B[0;32m~/Projects/COMO/main/como/rnaseq_gen.py:175\u001B[0m, in \u001B[0;36m_build_matrix_results\u001B[0;34m(matrix, gene_info, metadata_df, taxon)\u001B[0m\n\u001B[1;32m    167\u001B[0m \u001B[38;5;250m\u001B[39m\u001B[38;5;124;03m\"\"\"Read the counts matrix and returns the results.\u001B[39;00m\n\u001B[1;32m    168\u001B[0m \n\u001B[1;32m    169\u001B[0m \u001B[38;5;124;03m:param matrix: The gene counts matrix to process\u001B[39;00m\n\u001B[0;32m   (...)\u001B[0m\n\u001B[1;32m    172\u001B[0m \u001B[38;5;124;03m:return: A dataclass `ReadMatrixResults`\u001B[39;00m\n\u001B[1;32m    173\u001B[0m \u001B[38;5;124;03m\"\"\"\u001B[39;00m\n\u001B[1;32m    174\u001B[0m gene_info \u001B[38;5;241m=\u001B[39m gene_info_migrations(gene_info)\n\u001B[0;32m--> 175\u001B[0m conversion \u001B[38;5;241m=\u001B[39m \u001B[38;5;28;01mawait\u001B[39;00m ensembl_to_gene_id_and_symbol(ids\u001B[38;5;241m=\u001B[39m\u001B[43mmatrix\u001B[49m\u001B[43m[\u001B[49m\u001B[38;5;124;43m\"\u001B[39;49m\u001B[38;5;124;43mensembl_gene_id\u001B[39;49m\u001B[38;5;124;43m\"\u001B[39;49m\u001B[43m]\u001B[49m\u001B[38;5;241m.\u001B[39mtolist(), taxon\u001B[38;5;241m=\u001B[39mtaxon)\n\u001B[1;32m    176\u001B[0m matrix \u001B[38;5;241m=\u001B[39m matrix\u001B[38;5;241m.\u001B[39mmerge(conversion, on\u001B[38;5;241m=\u001B[39m\u001B[38;5;124m\"\u001B[39m\u001B[38;5;124mensembl_gene_id\u001B[39m\u001B[38;5;124m\"\u001B[39m, how\u001B[38;5;241m=\u001B[39m\u001B[38;5;124m\"\u001B[39m\u001B[38;5;124mleft\u001B[39m\u001B[38;5;124m\"\u001B[39m)\n\u001B[1;32m    178\u001B[0m \u001B[38;5;66;03m# Only include Entrez and Ensembl Gene IDs that are present in `gene_info`\u001B[39;00m\n",
      "File \u001B[0;32m~/Projects/COMO/.venv/lib/python3.10/site-packages/pandas/core/frame.py:4102\u001B[0m, in \u001B[0;36mDataFrame.__getitem__\u001B[0;34m(self, key)\u001B[0m\n\u001B[1;32m   4100\u001B[0m \u001B[38;5;28;01mif\u001B[39;00m \u001B[38;5;28mself\u001B[39m\u001B[38;5;241m.\u001B[39mcolumns\u001B[38;5;241m.\u001B[39mnlevels \u001B[38;5;241m>\u001B[39m \u001B[38;5;241m1\u001B[39m:\n\u001B[1;32m   4101\u001B[0m     \u001B[38;5;28;01mreturn\u001B[39;00m \u001B[38;5;28mself\u001B[39m\u001B[38;5;241m.\u001B[39m_getitem_multilevel(key)\n\u001B[0;32m-> 4102\u001B[0m indexer \u001B[38;5;241m=\u001B[39m \u001B[38;5;28;43mself\u001B[39;49m\u001B[38;5;241;43m.\u001B[39;49m\u001B[43mcolumns\u001B[49m\u001B[38;5;241;43m.\u001B[39;49m\u001B[43mget_loc\u001B[49m\u001B[43m(\u001B[49m\u001B[43mkey\u001B[49m\u001B[43m)\u001B[49m\n\u001B[1;32m   4103\u001B[0m \u001B[38;5;28;01mif\u001B[39;00m is_integer(indexer):\n\u001B[1;32m   4104\u001B[0m     indexer \u001B[38;5;241m=\u001B[39m [indexer]\n",
      "File \u001B[0;32m~/Projects/COMO/.venv/lib/python3.10/site-packages/pandas/core/indexes/base.py:3812\u001B[0m, in \u001B[0;36mIndex.get_loc\u001B[0;34m(self, key)\u001B[0m\n\u001B[1;32m   3807\u001B[0m     \u001B[38;5;28;01mif\u001B[39;00m \u001B[38;5;28misinstance\u001B[39m(casted_key, \u001B[38;5;28mslice\u001B[39m) \u001B[38;5;129;01mor\u001B[39;00m (\n\u001B[1;32m   3808\u001B[0m         \u001B[38;5;28misinstance\u001B[39m(casted_key, abc\u001B[38;5;241m.\u001B[39mIterable)\n\u001B[1;32m   3809\u001B[0m         \u001B[38;5;129;01mand\u001B[39;00m \u001B[38;5;28many\u001B[39m(\u001B[38;5;28misinstance\u001B[39m(x, \u001B[38;5;28mslice\u001B[39m) \u001B[38;5;28;01mfor\u001B[39;00m x \u001B[38;5;129;01min\u001B[39;00m casted_key)\n\u001B[1;32m   3810\u001B[0m     ):\n\u001B[1;32m   3811\u001B[0m         \u001B[38;5;28;01mraise\u001B[39;00m InvalidIndexError(key)\n\u001B[0;32m-> 3812\u001B[0m     \u001B[38;5;28;01mraise\u001B[39;00m \u001B[38;5;167;01mKeyError\u001B[39;00m(key) \u001B[38;5;28;01mfrom\u001B[39;00m \u001B[38;5;21;01merr\u001B[39;00m\n\u001B[1;32m   3813\u001B[0m \u001B[38;5;28;01mexcept\u001B[39;00m \u001B[38;5;167;01mTypeError\u001B[39;00m:\n\u001B[1;32m   3814\u001B[0m     \u001B[38;5;66;03m# If we have a listlike key, _check_indexing_error will raise\u001B[39;00m\n\u001B[1;32m   3815\u001B[0m     \u001B[38;5;66;03m#  InvalidIndexError. Otherwise we fall through and re-raise\u001B[39;00m\n\u001B[1;32m   3816\u001B[0m     \u001B[38;5;66;03m#  the TypeError.\u001B[39;00m\n\u001B[1;32m   3817\u001B[0m     \u001B[38;5;28mself\u001B[39m\u001B[38;5;241m.\u001B[39m_check_indexing_error(key)\n",
      "\u001B[0;31mKeyError\u001B[0m: 'ensembl_gene_id'"
     ]
    }
   ],
   "execution_count": 15
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Single-Cell RNA Sequencing Generation\n",
    "These parameters are identical to the ones listed for [total RNA sequencing](#Total-RNA-Sequencing-Generation), but they are listed again here for ease of reference\n",
    "\n",
    "### Parameters\n",
    "- `scrnaseq_config_file`: The configuration filename for total RNA. This file is found under the `data/config_sheets` folder\n",
    "- `rep_ratio`: The proportion of replicates before a gene is considered \"active\" in a study\n",
    "- `group_ratio`: The proportion of studies with expression required for a gene to be considered \"active\"\n",
    "- `rep_ratio_h`: The proportion of replicates that must express a gene before that gene is considered \"high-confidience\"\n",
    "- `group_ratio_h`: The proportion of studies that must express a gene before that gene is considered \"high-confidence\"\n",
    "- `technique`: The only option offered for single-cell RNA sequencing us \"umi\"\n",
    "- `quantile`: The cutoff Transcripts-Per-Million quantile for filtering\n",
    "- `min_zfpkm`: The cutoff for Counts-Per-Million filtering\n",
    "- `prep_method`: The library method used for preparation. Options are: `\"total\"`, `\"mRNA\"`, or `\"scrna\"`,\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from como.rnaseq_gen import FilteringTechnique, rnaseq_gen\n",
    "\n",
    "replicate_ratio = 0.75\n",
    "high_confidence_replicate_ratio = 1.0\n",
    "batch_ratio = 0.75\n",
    "high_confidence_batch_ratio = 1.0\n",
    "technique = FilteringTechnique.UMI\n",
    "cutoff = -3\n",
    "\n",
    "for context in context_names:\n",
    "    if context not in scrna_matrix_filepath:\n",
    "        continue\n",
    "    await rnaseq_gen(\n",
    "        context_name=context,\n",
    "        input_rnaseq_filepath=scrna_matrix_filepath[context],\n",
    "        input_gene_info_filepath=gene_info_filepath[context],\n",
    "        output_boolean_activity_filepath=scrna_matrix_filepath[context],\n",
    "        prep=RNAType.SCRNA,\n",
    "        taxon_id=taxon_id,\n",
    "        input_metadata_filepath=Path(\"./data/config_sheets/scrna_config.xlsx\"),\n",
    "        replicate_ratio=replicate_ratio,\n",
    "        high_replicate_ratio=high_confidence_replicate_ratio,\n",
    "        batch_ratio=batch_ratio,\n",
    "        high_batch_ratio=high_confidence_batch_ratio,\n",
    "        technique=technique,\n",
    "        cutoff=cutoff\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Proteomics Analysis\n",
    "The parameters here are mostly the same to total RNA and mRNA sequencing analysis, and are listed here for easier reference\n",
    "\n",
    "### Parameters\n",
    "- `proteomics_config_file`: The file path to the proteomics configuration file\n",
    "- `rep_ratio`: The ratio required before a gene is considered active in the replicate\n",
    "- `batch_ratio`: The ratio required before a gene is considered active in the study\n",
    "- `high_rep_ratio`: The ratio required before a gene is considered \"high-confidence\" in the replicate\n",
    "- `high_batch_ratio`: The ratio required before a gene is considered \"high-confidence\" in the study\n",
    "- `quantile`: The cutoff Transcripts-Per-Million quantile for filtering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from como.proteomics_gen import proteomics_gen\n",
    "\n",
    "for context in context_names:\n",
    "    await proteomics_gen(\n",
    "        context_name=context_names,\n",
    "        config_filepath=proteomics_metadata_filepath,\n",
    "        matrix_filepath=proteomics_matrix_filepath[context],\n",
    "        output_boolean_filepath=Path(f\"data/results/{context}/proteomics/{context}_proteomics_boolean_matrix.csv\"),\n",
    "        output_gaussian_png_filepath=Path(f\"data/results/{context}/proteomics/{context}_proteomics_gaussian.png\"),\n",
    "        output_gaussian_html_filepath=Path(f\"data/results/{context}/proteomics/{context}_proteomics_gaussian.html\"),\n",
    "        output_z_score_matrix_filepath=Path(f\"data/results/{context}/proteomics/{context}_zscore_matrix.csv\"),\n",
    "        input_entrez_map=Path(f\"data/results/{context}/proteomics/{context}_entrez_map.csv\"),\n",
    "        replicate_ratio=0.5,\n",
    "        batch_ratio=0.5,\n",
    "        high_confidence_replicate_ratio=0.7,\n",
    "        high_confidence_batch_ratio=0.7,\n",
    "        quantile=25,\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Cluster Sample Data (Optional)\n",
    "This step is used to cluster the samples based on their expression values. This can be used to determine which samples are more similar to each other. In a perfect world, one cluster would be created for each context type used. This is done using the `como/cluster_rnaseq.py` script.\n",
    "\n",
    "To see more about clustering, please visit the [Wikipedia article](https://en.wikipedia.org/wiki/Cluster_analysis)\n",
    "\n",
    "The parameters for this script are as follows:\n",
    "- `context_names`: The tissue/cell name of models that should be clustered. This was defined in the first code block, so it is not redefined here\n",
    "- `filt_technique`: The filtering technique to use; options are: `\"zfpkm\"`, `\"quantile\"`, or `\"cpm\"`\n",
    "- `cluster_algorithm`: The clustering algorithm to use. Options are: `\"mca\"` or `\"umap\"`\n",
    "- `label`: Should the samples be labeled in the plot? Options are: `\"True\"` or `\"False\"`\n",
    "- `min_dist`: The minimum distance for UMAP clustering. Must be between 0 and 1. Default value is 0.01\n",
    "- `replicate_ratio`: The ratio of active genes in replicates for a batch/study to be considered active. The default is 0.9\n",
    "- `batch_ratio`: The ratio of active genes in batches/studies for a context to be considered active. The default is 0.9\n",
    "- `min_count`: The ratio of active genes in a batch/study for a context to be considered active. The default is `\"default\"`\n",
    "- `quantile`: The ratio of active genes in a batch/study for a context to be considered active. The default is 0.5\n",
    "- `n_neighbors_rep`: N nearest neighbors for replicate clustering. The default is `\"default\"`, which is the total number of replicates\n",
    "- `n_neighbors_batch`: N nearest neighbors for batch clustering. The default is `\"default\"`, which is the total number of batches\n",
    "- `n_neighbors_context`: N nearest neighbors for context clustering. The default is `\"default\"`, which is the total number of contexts\n",
    "- `seed`: The random seed for clustering algorithm initialization. If not specified, `np.random.randint(0, 100000)` is used"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "filt_technique = \"zfpkm\"\n",
    "cluster_algorithm = \"umap\"\n",
    "label = True\n",
    "min_dist = 0.01\n",
    "replicate_ratio = 0.9\n",
    "batch_ratio = 0.9\n",
    "min_count = \"default\"\n",
    "quantile = 50\n",
    "n_neighbors_rep = \"default\"\n",
    "n_neighbors_batch = \"default\"\n",
    "n_neighbors_context = \"default\"\n",
    "seed = -1\n",
    "\n",
    "# fmt: off\n",
    "cmd = \" \".join(\n",
    "    [\n",
    "        \"python3\", \"como/cluster_rnaseq.py\",\n",
    "        \"--context-names\", context_names,\n",
    "        \"--filt-technique\", filt_technique,\n",
    "        \"--cluster-algorithm\", cluster_algorithm,\n",
    "        \"--label\", label,\n",
    "        \"--min-dist\", str(min_dist),\n",
    "        \"--replicate-ratio\", str(replicate_ratio),\n",
    "        \"--batch-ratio\", str(batch_ratio),\n",
    "        \"--n-neighbors-rep\", str(n_neighbors_rep),\n",
    "        \"--n-neighbors-batch\", str(n_neighbors_batch),\n",
    "        \"--n-neighbors-context\", str(n_neighbors_context),\n",
    "        \"--min-count\", str(min_count),\n",
    "        \"--quantile\", str(quantile),\n",
    "        \"--seed\", str(seed),\n",
    "    ]\n",
    ")\n",
    "# fmt: on\n",
    "\n",
    "!{cmd}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Merge Expression from Different Data Sources\n",
    "\n",
    "Thus far, active genes have been determined for at least one data source. If multiple data sources are being used, we can merge the active genes from these sources to make a list of active genes that is more comprehensive (or strict!) than any data source on its own.\n",
    "\n",
    "`como/merge_xomics.py` takes each data source discussed so far as an argument. The other arguments to consider are:\n",
    "- `--expression-requirement`: The number of data sources with expression required for a gene to be considered active, if the gene is not \"high-confidence\" for any source. (default: total number of input sources provided)\n",
    "- `--requirement-adjust`: This is used to adjust the expression requirement argument in the event that tissues have a different number of provided data sources. This does nothing if there is only one tissue type in the configuration files.\n",
    "    - `\"progressive\"`: The expression requirement applies to tissue(s) with the lowest number of data sources. Tissues with more than this value will require its genes to be expressed in 1 additional source before it is \"active\" in the model\n",
    "    - `\"regressive\"` (default): The expression requirement applies to the tissue(s) with the largest number of data sources. Tissues with less than this value will require its genes to be expressed in 1 fewer sources before the gene is considered \"active\" in the model.\n",
    "    - `\"flat\"`: The expression requirement is used regardless of differences in the number of data sources provided for different tissues\n",
    "\n",
    "- `--no-hc`: This flag should be set to prevent high-confidence genes from overriding the expression requirement set.\n",
    "    - If this flag is not used, any gene that was determined to be \"high-confidence\" in any input source will cause the gene to be active in the final model, regardless of agreement with other sources\n",
    "- `--no-na-adjustment`: This flag should be used to prevent genes that are not present in one data source, but are present in others, from subtracting one from the expression requirement.\n",
    "    - If this flag is not used, any time a gene is \"NA\" in a source, meaning it was not tested for in the library of that data sources but <ins>was</ins> tested in the library of another source, it will subtract one from the expression requirement.\n",
    "\n",
    "The adjusted expression requirement will never resolve to be less than one or greater than the number of data sources for a given tissue\n",
    "\n",
    "### Parameters\n",
    "The three parameters listed here were used in RNA Sequencing generation, and should not need to be defined. If you did **not** use one of these, simply un-comment it from the command below by placing a \"`#`\" at the beginning of the appropriate lines\n",
    "- `trnaseq_config_file`: The file name used in the [total RNA Sequencing](#Total-RNA-Sequencing-Generation) section of the notebook\n",
    "- `mrnaseq_config_file`: The file name used in the [mRNA Sequencing](#mRNA-Sequencing-Generation) section of the notebook\n",
    "- `proteomics_config_file`: The file name used in the [proteomics generation](#Proteomics-Analysis) section of the notebook\n",
    "\n",
    "The following parameters have not been used in a previous section of the notebook, so they are defined in the below code block\n",
    "- `expression_requirement`: This is the number of sources a gene must be active in for it to be considered active\n",
    "- `requirement_adjust`: The technique to adjust expression requirement based on differences in number of provided data source types\n",
    "- `total_rna_weight`: Total RNA-seq weight for merging zFPKM distribution\n",
    "- `mrna_weight`: mRNA weight for merging zFPKM distribution\n",
    "- `single_cell_weight`: Single-cell weight for merging zFPKM distribution\n",
    "- `proteomics_weight`: Proteomic weight for merging zFPKM distribution\n",
    "\n",
    "Each of the \"weights\" (`total_rna_weight`, `mrna_weight`, etc.) are used to place a significance on each method. Becuase there are many steps in the Dogma from transcription to translation, the gene expression as seen by total RNA or mRNA sequencing may not be representative of the gene's protein expression, and this its metabolic impact. Because of this, you are able to weight each source more (or less) than another."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-12-16T19:54:26.540397Z",
     "start_time": "2024-12-16T19:54:26.515554Z"
    }
   },
   "source": [
    "from como.merge_xomics import AdjustmentMethod, merge_xomics\n",
    "\n",
    "total_rna_weight = 6\n",
    "mrna_weight = 6\n",
    "single_cell_weight = 6\n",
    "proteomics_weight = 10\n",
    "minimum_source_expression = 1\n",
    "expression_requirement = 1\n",
    "requirement_adjustment_method = AdjustmentMethod.REGRESSIVE\n",
    "force_activate_high_confidence = False\n",
    "adjust_for_na_sources = False\n",
    "merge_zfpkm_distrubution = True\n",
    "keep_transcriptomics_score = True\n",
    "\n",
    "\n",
    "mrna_batches = {\"naiveB\": [\n",
    "    \"naiveB_S1R1\", \"naiveB_S1R2\", \"naiveB_S1R3\", \"naiveB_S1R4\",\n",
    "    \"naiveB_S2R1\", \"naiveB_S2R2\", \"naiveB_S2R3\", \"naiveB_S2R4\",\n",
    "]}\n",
    "trna_batches = {\"naiveB\": [\"naiveB_S3R1\", \"naiveB_S3R2\", \"naiveB_S3R3\"]}\n",
    "\n",
    "\n",
    "for context in context_names:\n",
    "    await merge_xomics(\n",
    "        context_name=context,\n",
    "        trna_matrix_or_filepath=trna_matrix_filepath[context],\n",
    "        mrna_matrix_or_filepath=mrna_matrix_filepath[context],\n",
    "        scrna_matrix_or_filepath=None,  # scrna_matrix_filepath[context],\n",
    "        proteomic_matrix_or_filepath=None,  # proteomics_matrix_filepath[context],\n",
    "        trna_batches=trna_batches,\n",
    "        mrna_batches=mrna_batches,\n",
    "        scrna_batches=None,\n",
    "        proteomic_batches=None,\n",
    "        trna_weight=total_rna_weight,\n",
    "        mrna_weight=mrna_weight,\n",
    "        scrna_weight=single_cell_weight,\n",
    "        proteomic_weight=proteomics_weight,\n",
    "        minimum_source_expression=minimum_source_expression,\n",
    "        expression_requirement=expression_requirement,\n",
    "        adjust_method=requirement_adjustment_method,\n",
    "        force_activate_high_confidence=force_activate_high_confidence,\n",
    "        adjust_for_na=adjust_for_na_sources,\n",
    "        merge_zfpkm_distribution=merge_zfpkm_distrubution,\n",
    "        keep_transcriptomics_score=keep_transcriptomics_score,\n",
    "        output_merge_activity_filepath=Path(f\"data/results/{context}/ActiveGenes_{context}_Merged.csv\"),\n",
    "        output_transcriptomic_details_filepath=Path(f\"data/results/{context}/TranscriptomicDetails_{context}.csv\"),\n",
    "        output_trna_activity_filepath=Path(f\"data/results/{context}/total-rna/trna_activity_{context}.csv\"),\n",
    "        output_mrna_activity_filepath=Path(f\"data/results/{context}/mrna/mrna_activity_{context}.csv\"),\n",
    "        output_scrna_activity_filepath=Path(f\"data/results/{context}/scrna/scrna_activity_{context}.csv\"),\n",
    "        output_proteomic_activity_filepath=Path(f\"data/results/{context}/proteomics/proteomic_activity_{context}.csv\"),\n",
    "        output_final_model_scores_filepath=Path(f\"data/results/{context}/model_scores_{context}.csv\"),\n",
    "        output_figure_dirpath=Path(f\"data/results/{context}/figures\")\n",
    "    )"
   ],
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "data/results/naiveB/total-rna/trna_activity_naiveB.csv\n"
     ]
    },
    {
     "ename": "FileExistsError",
     "evalue": "[Errno 17] File exists: 'data/results/naiveB/total-rna'",
     "output_type": "error",
     "traceback": [
      "\u001B[0;31m---------------------------------------------------------------------------\u001B[0m",
      "\u001B[0;31mFileExistsError\u001B[0m                           Traceback (most recent call last)",
      "Cell \u001B[0;32mIn[6], line 24\u001B[0m\n\u001B[1;32m     20\u001B[0m trna_batches \u001B[38;5;241m=\u001B[39m {\u001B[38;5;124m\"\u001B[39m\u001B[38;5;124mnaiveB\u001B[39m\u001B[38;5;124m\"\u001B[39m: [\u001B[38;5;124m\"\u001B[39m\u001B[38;5;124mnaiveB_S3R1\u001B[39m\u001B[38;5;124m\"\u001B[39m, \u001B[38;5;124m\"\u001B[39m\u001B[38;5;124mnaiveB_S3R2\u001B[39m\u001B[38;5;124m\"\u001B[39m, \u001B[38;5;124m\"\u001B[39m\u001B[38;5;124mnaiveB_S3R3\u001B[39m\u001B[38;5;124m\"\u001B[39m]}\n\u001B[1;32m     23\u001B[0m \u001B[38;5;28;01mfor\u001B[39;00m context \u001B[38;5;129;01min\u001B[39;00m context_names:\n\u001B[0;32m---> 24\u001B[0m     \u001B[38;5;28;01mawait\u001B[39;00m merge_xomics(\n\u001B[1;32m     25\u001B[0m         context_name\u001B[38;5;241m=\u001B[39mcontext_names,\n\u001B[1;32m     26\u001B[0m         trna_matrix_or_filepath\u001B[38;5;241m=\u001B[39mtrna_matrix_filepath[context],\n\u001B[1;32m     27\u001B[0m         mrna_matrix_or_filepath\u001B[38;5;241m=\u001B[39mmrna_matrix_filepath[context],\n\u001B[1;32m     28\u001B[0m         scrna_matrix_or_filepath\u001B[38;5;241m=\u001B[39m\u001B[38;5;28;01mNone\u001B[39;00m,  \u001B[38;5;66;03m# scrna_matrix_filepath[context],\u001B[39;00m\n\u001B[1;32m     29\u001B[0m         proteomic_matrix_or_filepath\u001B[38;5;241m=\u001B[39m\u001B[38;5;28;01mNone\u001B[39;00m,  \u001B[38;5;66;03m# proteomics_matrix_filepath[context],\u001B[39;00m\n\u001B[1;32m     30\u001B[0m         trna_batches\u001B[38;5;241m=\u001B[39mtrna_batches,\n\u001B[1;32m     31\u001B[0m         mrna_batches\u001B[38;5;241m=\u001B[39mmrna_batches,\n\u001B[1;32m     32\u001B[0m         scrna_batches\u001B[38;5;241m=\u001B[39m\u001B[38;5;28;01mNone\u001B[39;00m,\n\u001B[1;32m     33\u001B[0m         proteomic_batches\u001B[38;5;241m=\u001B[39m\u001B[38;5;28;01mNone\u001B[39;00m,\n\u001B[1;32m     34\u001B[0m         trna_weight\u001B[38;5;241m=\u001B[39mtotal_rna_weight,\n\u001B[1;32m     35\u001B[0m         mrna_weight\u001B[38;5;241m=\u001B[39mmrna_weight,\n\u001B[1;32m     36\u001B[0m         scrna_weight\u001B[38;5;241m=\u001B[39msingle_cell_weight,\n\u001B[1;32m     37\u001B[0m         proteomic_weight\u001B[38;5;241m=\u001B[39mproteomics_weight,\n\u001B[1;32m     38\u001B[0m         minimum_source_expression\u001B[38;5;241m=\u001B[39mminimum_source_expression,\n\u001B[1;32m     39\u001B[0m         expression_requirement\u001B[38;5;241m=\u001B[39mexpression_requirement,\n\u001B[1;32m     40\u001B[0m         adjust_method\u001B[38;5;241m=\u001B[39mrequirement_adjustment_method,\n\u001B[1;32m     41\u001B[0m         force_activate_high_confidence\u001B[38;5;241m=\u001B[39mforce_activate_high_confidence,\n\u001B[1;32m     42\u001B[0m         adjust_for_na\u001B[38;5;241m=\u001B[39madjust_for_na_sources,\n\u001B[1;32m     43\u001B[0m         merge_zfpkm_distribution\u001B[38;5;241m=\u001B[39mmerge_zfpkm_distrubution,\n\u001B[1;32m     44\u001B[0m         keep_transcriptomics_score\u001B[38;5;241m=\u001B[39mkeep_transcriptomics_score,\n\u001B[1;32m     45\u001B[0m         output_merge_activity_filepath\u001B[38;5;241m=\u001B[39mPath(\u001B[38;5;124mf\u001B[39m\u001B[38;5;124m\"\u001B[39m\u001B[38;5;124mdata/results/\u001B[39m\u001B[38;5;132;01m{\u001B[39;00mcontext\u001B[38;5;132;01m}\u001B[39;00m\u001B[38;5;124m/ActiveGenes_\u001B[39m\u001B[38;5;132;01m{\u001B[39;00mcontext\u001B[38;5;132;01m}\u001B[39;00m\u001B[38;5;124m_Merged.csv\u001B[39m\u001B[38;5;124m\"\u001B[39m),\n\u001B[1;32m     46\u001B[0m         output_transcriptomic_details_filepath\u001B[38;5;241m=\u001B[39mPath(\u001B[38;5;124mf\u001B[39m\u001B[38;5;124m\"\u001B[39m\u001B[38;5;124mdata/results/\u001B[39m\u001B[38;5;132;01m{\u001B[39;00mcontext\u001B[38;5;132;01m}\u001B[39;00m\u001B[38;5;124m/TranscriptomicDetails_\u001B[39m\u001B[38;5;132;01m{\u001B[39;00mcontext\u001B[38;5;132;01m}\u001B[39;00m\u001B[38;5;124m.csv\u001B[39m\u001B[38;5;124m\"\u001B[39m),\n\u001B[1;32m     47\u001B[0m         output_trna_activity_filepath\u001B[38;5;241m=\u001B[39mPath(\u001B[38;5;124mf\u001B[39m\u001B[38;5;124m\"\u001B[39m\u001B[38;5;124mdata/results/\u001B[39m\u001B[38;5;132;01m{\u001B[39;00mcontext\u001B[38;5;132;01m}\u001B[39;00m\u001B[38;5;124m/total-rna/trna_activity_\u001B[39m\u001B[38;5;132;01m{\u001B[39;00mcontext\u001B[38;5;132;01m}\u001B[39;00m\u001B[38;5;124m.csv\u001B[39m\u001B[38;5;124m\"\u001B[39m),\n\u001B[1;32m     48\u001B[0m         output_mrna_activity_filepath\u001B[38;5;241m=\u001B[39mPath(\u001B[38;5;124mf\u001B[39m\u001B[38;5;124m\"\u001B[39m\u001B[38;5;124mdata/results/\u001B[39m\u001B[38;5;132;01m{\u001B[39;00mcontext\u001B[38;5;132;01m}\u001B[39;00m\u001B[38;5;124m/mrna/mrna_activity_\u001B[39m\u001B[38;5;132;01m{\u001B[39;00mcontext\u001B[38;5;132;01m}\u001B[39;00m\u001B[38;5;124m.csv\u001B[39m\u001B[38;5;124m\"\u001B[39m),\n\u001B[1;32m     49\u001B[0m         output_scrna_activity_filepath\u001B[38;5;241m=\u001B[39mPath(\u001B[38;5;124mf\u001B[39m\u001B[38;5;124m\"\u001B[39m\u001B[38;5;124mdata/results/\u001B[39m\u001B[38;5;132;01m{\u001B[39;00mcontext\u001B[38;5;132;01m}\u001B[39;00m\u001B[38;5;124m/scrna/scrna_activity_\u001B[39m\u001B[38;5;132;01m{\u001B[39;00mcontext\u001B[38;5;132;01m}\u001B[39;00m\u001B[38;5;124m.csv\u001B[39m\u001B[38;5;124m\"\u001B[39m),\n\u001B[1;32m     50\u001B[0m         output_proteomic_activity_filepath\u001B[38;5;241m=\u001B[39mPath(\u001B[38;5;124mf\u001B[39m\u001B[38;5;124m\"\u001B[39m\u001B[38;5;124mdata/results/\u001B[39m\u001B[38;5;132;01m{\u001B[39;00mcontext\u001B[38;5;132;01m}\u001B[39;00m\u001B[38;5;124m/proteomics/proteomic_activity_\u001B[39m\u001B[38;5;132;01m{\u001B[39;00mcontext\u001B[38;5;132;01m}\u001B[39;00m\u001B[38;5;124m.csv\u001B[39m\u001B[38;5;124m\"\u001B[39m),\n\u001B[1;32m     51\u001B[0m         output_final_model_scores_filepath\u001B[38;5;241m=\u001B[39mPath(\u001B[38;5;124mf\u001B[39m\u001B[38;5;124m\"\u001B[39m\u001B[38;5;124mdata/results/\u001B[39m\u001B[38;5;132;01m{\u001B[39;00mcontext\u001B[38;5;132;01m}\u001B[39;00m\u001B[38;5;124m/model_scores_\u001B[39m\u001B[38;5;132;01m{\u001B[39;00mcontext\u001B[38;5;132;01m}\u001B[39;00m\u001B[38;5;124m.csv\u001B[39m\u001B[38;5;124m\"\u001B[39m),\n\u001B[1;32m     52\u001B[0m         output_figure_dirpath\u001B[38;5;241m=\u001B[39mPath(\u001B[38;5;124mf\u001B[39m\u001B[38;5;124m\"\u001B[39m\u001B[38;5;124mdata/results/\u001B[39m\u001B[38;5;132;01m{\u001B[39;00mcontext\u001B[38;5;132;01m}\u001B[39;00m\u001B[38;5;124m/figures\u001B[39m\u001B[38;5;124m\"\u001B[39m)\n\u001B[1;32m     53\u001B[0m     )\n",
      "File \u001B[0;32m~/Projects/COMO/main/como/merge_xomics.py:678\u001B[0m, in \u001B[0;36mmerge_xomics\u001B[0;34m(context_name, trna_matrix_or_filepath, mrna_matrix_or_filepath, scrna_matrix_or_filepath, proteomic_matrix_or_filepath, trna_batches, mrna_batches, scrna_batches, proteomic_batches, output_merge_activity_filepath, output_transcriptomic_details_filepath, output_trna_activity_filepath, output_mrna_activity_filepath, output_scrna_activity_filepath, output_proteomic_activity_filepath, output_final_model_scores_filepath, output_figure_dirpath, trna_weight, mrna_weight, scrna_weight, proteomic_weight, minimum_source_expression, expression_requirement, adjust_method, force_activate_high_confidence, adjust_for_na, merge_zfpkm_distribution, keep_transcriptomics_score, weighted_z_floor, weighted_z_ceiling)\u001B[0m\n\u001B[1;32m    676\u001B[0m     output_transcriptomic_details_filepath\u001B[38;5;241m.\u001B[39mparent\u001B[38;5;241m.\u001B[39mmkdir(parents\u001B[38;5;241m=\u001B[39m\u001B[38;5;28;01mTrue\u001B[39;00m)\n\u001B[1;32m    677\u001B[0m \u001B[38;5;28;01mif\u001B[39;00m output_trna_activity_filepath \u001B[38;5;129;01mand\u001B[39;00m \u001B[38;5;129;01mnot\u001B[39;00m output_trna_activity_filepath\u001B[38;5;241m.\u001B[39mexists():\n\u001B[0;32m--> 678\u001B[0m     \u001B[43moutput_trna_activity_filepath\u001B[49m\u001B[38;5;241;43m.\u001B[39;49m\u001B[43mparent\u001B[49m\u001B[38;5;241;43m.\u001B[39;49m\u001B[43mmkdir\u001B[49m\u001B[43m(\u001B[49m\u001B[43mparents\u001B[49m\u001B[38;5;241;43m=\u001B[39;49m\u001B[38;5;28;43;01mTrue\u001B[39;49;00m\u001B[43m)\u001B[49m\n\u001B[1;32m    679\u001B[0m \u001B[38;5;28;01mif\u001B[39;00m output_mrna_activity_filepath \u001B[38;5;129;01mand\u001B[39;00m \u001B[38;5;129;01mnot\u001B[39;00m output_mrna_activity_filepath\u001B[38;5;241m.\u001B[39mexists():\n\u001B[1;32m    680\u001B[0m     output_mrna_activity_filepath\u001B[38;5;241m.\u001B[39mparent\u001B[38;5;241m.\u001B[39mmkdir(parents\u001B[38;5;241m=\u001B[39m\u001B[38;5;28;01mTrue\u001B[39;00m)\n",
      "File \u001B[0;32m~/.local/share/uv/python/cpython-3.10.15-macos-aarch64-none/lib/python3.10/pathlib.py:1175\u001B[0m, in \u001B[0;36mPath.mkdir\u001B[0;34m(self, mode, parents, exist_ok)\u001B[0m\n\u001B[1;32m   1171\u001B[0m \u001B[38;5;250m\u001B[39m\u001B[38;5;124;03m\"\"\"\u001B[39;00m\n\u001B[1;32m   1172\u001B[0m \u001B[38;5;124;03mCreate a new directory at this given path.\u001B[39;00m\n\u001B[1;32m   1173\u001B[0m \u001B[38;5;124;03m\"\"\"\u001B[39;00m\n\u001B[1;32m   1174\u001B[0m \u001B[38;5;28;01mtry\u001B[39;00m:\n\u001B[0;32m-> 1175\u001B[0m     \u001B[38;5;28;43mself\u001B[39;49m\u001B[38;5;241;43m.\u001B[39;49m\u001B[43m_accessor\u001B[49m\u001B[38;5;241;43m.\u001B[39;49m\u001B[43mmkdir\u001B[49m\u001B[43m(\u001B[49m\u001B[38;5;28;43mself\u001B[39;49m\u001B[43m,\u001B[49m\u001B[43m \u001B[49m\u001B[43mmode\u001B[49m\u001B[43m)\u001B[49m\n\u001B[1;32m   1176\u001B[0m \u001B[38;5;28;01mexcept\u001B[39;00m \u001B[38;5;167;01mFileNotFoundError\u001B[39;00m:\n\u001B[1;32m   1177\u001B[0m     \u001B[38;5;28;01mif\u001B[39;00m \u001B[38;5;129;01mnot\u001B[39;00m parents \u001B[38;5;129;01mor\u001B[39;00m \u001B[38;5;28mself\u001B[39m\u001B[38;5;241m.\u001B[39mparent \u001B[38;5;241m==\u001B[39m \u001B[38;5;28mself\u001B[39m:\n",
      "\u001B[0;31mFileExistsError\u001B[0m: [Errno 17] File exists: 'data/results/naiveB/total-rna'"
     ]
    }
   ],
   "execution_count": 6
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Step 2: Create Tissue/Cell-Type Specific Models\n",
    "\n",
    "## Boundary Reactions\n",
    "To create a metabolic model, the following information about each metabolite or reaction involved is required:\n",
    "- **Reaction Type**\n",
    "    - Exchange\n",
    "    - Demand\n",
    "    - Sink\n",
    "- **Metabolic/Reaction Abbreviation**\n",
    "    - You can use the [Virutal Metabolic Human](https://www.vmh.life/#home) to look up your metabolite and reaction abbreviations\n",
    "- **Compartments**\n",
    "    - Cytosol\n",
    "    - Extracellular\n",
    "    - Golgi Apparatus\n",
    "    - Internal Membranes\n",
    "    - Lysosome\n",
    "    - Mitochondria\n",
    "    - Nucleus\n",
    "    - Endoplasmic Reticulum\n",
    "    - Unknown\n",
    "- **Minimum Reaction Rate**\n",
    "- **Maximum Reaction Rate**\n",
    "\n",
    "\n",
    "*Below is an example of a properly formatted table of metabolic and reaction information*\n",
    "\n",
    "| Reaction | Abbreviation |    Compartment     | Minimum Reaction Rate | Maximum Reaction Rate |\n",
    "|:--------:|:------------:|:------------------:|:---------------------:|:---------------------:|\n",
    "| Exchange |    glc_D     |   Extracellular    |         -100          |         1000          |\n",
    "|  Demand  |  15HPETATP   |      Cytosol       |          -1           |         1000          |\n",
    "|   Sink   |    met_L     | Internal Membranes |         -1000         |           1           |\n",
    "\n",
    "\n",
    "These reactions should be placed into a CSV file; a template can be found at `data/boundary_rxns/default_force_rxns.csv`. Append your reactions to this file, and remove any that are not required. COMO will load this file in during model creation\n",
    "\n",
    "## Force Reactions\n",
    "Force reactions are reactions that should **always** be included in the model, no matter their flux value in the metabolic data provided. In contrast to the boundary reaction list, this is simply a list of reaction names that should be \"forced\" through the model. Append your force reactions to the `data/force_rxns/default_force_rxns.csv` file, and remove any that are not required. COMO will load this file during model creation\n",
    "\n",
    "*Below is an example of a properly formatted table of force reactions*\n",
    "\n",
    "| Reaction |\n",
    "|:--------:|\n",
    "|  glc_D   |\n",
    "|  met_L   |\n",
    "\n",
    "## Adding Reference Models\n",
    "This Jupyter notebook uses Recon3D's [Virtual Metabolic Human](https://www.vmh.life/) as a base to map reactions onto, and is included with the Jupyter notebook. If you would like to include other reference models, simply upload them to the `data` folder, and set the name of the `general_model_file` below to the name of your reference model.\n",
    "\n",
    "## Parameters\n",
    "The following is a list of parameters and their function in this section of the pipeline\n",
    "- `low_thres`: If you are using the `IMAT` reconstruction algorithm, gene expression above this value will be placed in the \"mid-expression\" bin\n",
    "- `high_thres`: If you are using the `IMAT` reconstruction algorithm, gene expression above this value will be placed in the \"high-expression\" bin\n",
    "- `output_filetypes`: These are the file types you would like to save your model as. It should be one (or multiple) of the following: `\"xml\"`, `\"mat\"`, `\"json\"`\n",
    "- `objective_dict`: This is an objective the model should be solved for. Popular options are `\"biomass_reaction\"` or `\"biomass_maintenance\"`\n",
    "- `general_model_file`: This is the reference model file to load\n",
    "- `recon_algorithm`: The troppo reconstruction algorithm to use. This should be one of the following: `\"FastCORE\"`, `\"CORDA\"`, `\"GIMME\"`, `\"tINIT\"`, `\"IMAT\"`\n",
    "- `solver`: The solver to use for optimizing the model. Options are: `\"GUROBI\"` or `\"GLPK\"`\n",
    "- `boundary_reactions_filename`: The filename of boundary reactions that should be used\n",
    "- `force_reactions_filename`: The filename of the force reactions to be used. Force reactions will (as the name implies) force the optimizer to use these reactions, **no matter their expression**\n",
    "- `exclude_reactions_filename`: The filename of reactions to exclude from the model, no matter their expression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import json\n",
    "from pathlib import Path\n",
    "\n",
    "from como.project import Config\n",
    "\n",
    "# Set your objectives before running!\n",
    "objective_dict = {\"naiveB\": \"biomass_maintenance\", \"smB\": \"biomass_maintenance\"}\n",
    "# -----------------\n",
    "\n",
    "low_threshold = -5\n",
    "high_threshold = -3\n",
    "output_filetypes = \"xml mat json\"\n",
    "general_model_file = \"GeneralModelUpdatedV2.mat\"\n",
    "recon_algorithms = [\"IMAT\"]\n",
    "solver = \"GUROBI\"\n",
    "\n",
    "config = Config()\n",
    "\n",
    "# Load the output of step 1, which is a dictionary that specifies the merged list of active Gene IDs for each tissue\n",
    "step1_results_file = config.data_dir / \"results\" \"step1_results_files.json\"\n",
    "with step1_results_file.open(\"r\") as json_file:\n",
    "    context_gene_exp = json.load(json_file)\n",
    "\n",
    "for recon_algorithm in recon_algorithms:\n",
    "    for context in context_gene_exp:\n",
    "        objective = objective_dict[context]\n",
    "\n",
    "        if recon_algorithm.upper() in [\"IMAT\", \"TINIT\"]:\n",
    "            active_genes_filepath = config.data_dir / \"results\" / context / f\"model_scores_{context}.csv\"\n",
    "        else:\n",
    "            gene_expression_file = context_gene_exp[context]\n",
    "            active_genes_filename = Path(gene_expression_file).name\n",
    "            active_genes_filepath = config.data_dir / \"results\" / context / active_genes_filename\n",
    "\n",
    "        general_model_filepath = config.data_dir / \"GeneralModelUpdatedV2.mat\"\n",
    "        boundary_reactions_filepath = config.data_dir / \"boundary_rxns\" / f\"{context}_boundary_rxns.csv\"\n",
    "        force_reactions_filepath = config.data_dir / \"force_rxns\" / f\"{context}_force_rxns.csv\"\n",
    "        exclude_reactions_filepath = config.data_dir / \"exclude_rxns\" / f\"{context}_exclude_rxns.csv\"\n",
    "\n",
    "        # fmt: off\n",
    "        cmd = \" \".join(\n",
    "            [\n",
    "                \"python3\", \"como/create_context_specific_model.py\",\n",
    "                \"--context\", context,\n",
    "                \"--reference-model-filepath\", general_model_filepath,\n",
    "                \"--active-genes-filepath\", active_genes_filepath,\n",
    "                \"--objective\", objective,\n",
    "                \"--boundary-reactions-filepath\", boundary_reactions_filepath,\n",
    "                # \"--exclude-reactions-filepath\", exclude_reactions_filepath,\n",
    "                \"--force-reactions-filepath\", force_reactions_filepath,\n",
    "                \"--algorithm\", recon_algorithm,\n",
    "                \"--low-threshold\", str(low_threshold),\n",
    "                \"--high-threshold\", str(high_threshold),\n",
    "                \"--solver\", solver,\n",
    "                \"--output-filetypes\", output_filetypes,\n",
    "            ]\n",
    "        )\n",
    "        # fmt: on\n",
    "        !{cmd}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Generate MEMOTE Reports\n",
    "> NOTE: This step is entirely optional\n",
    "\n",
    "MEMOTE is an open-source tool to automate the testing and reporting of metabolic models. This report is a detailed summary of the tests performed by MEMOTE on a given metabolic model (i.e., the one you just generated), along with the results and recommendations for improving the model. In order to create these reports, a metabolic \"map\" is required. Several of these are included in COMO, found under `data/maps/RECON1`. If you would like to add your own maps, they can be included in multiple places:\n",
    "1. If you have mapped a `local_files` directory to the container before starting, you can simply copy-and-paste them into the `local_files/maps` directory using the file browser of your computer. This is the most robust solution because the files will not be deleted by the container after it stops, or if it is updated in the future\n",
    "2. You can upload them to the Jupyter notebook under the `data/maps` directory. The code block below will search for any `.json` files that are not already included in the `map_dict` dictionary\n",
    "\n",
    "config.data_dir,\n",
    "                \"results\",\n",
    "                context,\n",
    "                \"figures\",\n",
    "                f\"{key}_map_{context}_{algorithm}.html\"\n",
    "\n",
    "The resulting MEMOTE reports will be saved to `data/results/exampleTissue/figures/mapName_map_exampleTissue_ALGORITHM.html`.\n",
    "\n",
    "- `mapName`: This is the name of the map file. In the `map_dict` dictionary below, this value would be `trypto`, `retinol`, etc.\n",
    "- `exampleTissue`: This is the name of the tissue context\n",
    "- `ALGORITHM`: This is the algorithm (`recon_algorithm`) used in the above model creation step\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cobra\n",
    "from como.project import Config\n",
    "from escher import Builder\n",
    "from loguru import logger\n",
    "\n",
    "config = Config()\n",
    "\n",
    "user_map_dir = config.data_dir / \"local_files\" / \"maps\"\n",
    "map_dict = {\n",
    "    \"trypto\": f\"{config.data_dir}/maps/RECON1/RECON1.tryptophan_metabolism.json\",\n",
    "    # \"lipid\": f\"{config.data_dir}/maps/RECON1/RECON1.\",  # Not present in COMO by default yet\n",
    "    \"retinol\": f\"{config.data_dir}/maps/RECON1/RECON1.inositol_retinol_metabolism.json\",\n",
    "    \"glyco\": f\"{config.data_dir}/maps/RECON1/RECON1.glycolysis_TCA_PPP.json\",\n",
    "    \"combined\": f\"{config.data_dir}/maps/RECON1/RECON1.combined.json\",\n",
    "    \"carbo\": f\"{config.data_dir}/maps/RECON1/RECON1.carbohydrate_metabolism.json\",\n",
    "    \"amino\": f\"{config.data_dir}/maps/RECON1/RECON1.amino_acid_partial_metabolism.json\",\n",
    "}\n",
    "\n",
    "# Collect files from user-input json maps\n",
    "index = 1\n",
    "for file in user_map_dir.glob(\"**/*.json\"):\n",
    "    map_dict[file.stem] = file\n",
    "    index += 1\n",
    "\n",
    "# Collect any additional maps under the `{config.data_dir}/maps/` directory\n",
    "for file in (config.data_dir / \"maps\").glob(\"**/*.json\"):\n",
    "    if file not in map_dict.values():\n",
    "        map_dict[file.stem] = file\n",
    "\n",
    "for recon_algorithm in recon_algorithms:\n",
    "    for context in context_gene_exp:\n",
    "        model_json = config.data_dir / \"results\" / context / f\"{context}_SpecificModel_{recon_algorithm}.json\"\n",
    "\n",
    "        logger.info(f\"Loading '{context}', this may take some time...\")\n",
    "        model = cobra.io.load_json_model(model_json)\n",
    "        for key in map_dict:\n",
    "            logger.info(f\"Running with: {key}\")\n",
    "            builder = Builder(map_json=str(map_dict[key]))\n",
    "            builder.model = model\n",
    "            solution = cobra.flux_analysis.pfba(model)\n",
    "            builder.reaction_data = solution.fluxes\n",
    "            builder.reaction_scale = [\n",
    "                {\"type\": \"min\", \"color\": \"#ff3300\", \"size\": 12},\n",
    "                {\"type\": \"q1\", \"color\": \"#ffc61a\", \"size\": 14},\n",
    "                {\"type\": \"median\", \"color\": \"#ffe700\", \"size\": 16},\n",
    "                {\"type\": \"q3\", \"color\": \"#4ffd3c\", \"size\": 18},\n",
    "                {\"type\": \"max\", \"color\": \"#3399ff\", \"size\": 20},\n",
    "            ]\n",
    "            builder.reaction_no_data_color = \"#8e8e8e\"\n",
    "\n",
    "            builder.save_html(\n",
    "                config.data_dir / \"results\" / context / \"figures\" / f\"{key}_map_{context}_{recon_algorithm}.html\"\n",
    "            )\n",
    "\n",
    "        out_dir = config.data_dir / \"results\" / context\n",
    "        report_file = out_dir / f\"memote_report_{context}_{recon_algorithm}.html\"\n",
    "        model_file = out_dir / f\"{context}_SpecificModel_{recon_algorithm}.xml\"\n",
    "        log_dir = out_dir / \"memote\"\n",
    "        log_file = log_dir / f\"{context}_{recon_algorithm}_memote.log\"\n",
    "\n",
    "        if not log_dir.exists():\n",
    "            log_dir.mkdir(parents=True, exist_ok=True)\n",
    "\n",
    "        cmd = \" \".join(\n",
    "            [\"memote\", \"report\", \"snapshot\", \"--filename\", f\"{report_file}\", f\"{model_file}\", \">\", f\"{log_file}\"]\n",
    "        )\n",
    "\n",
    "        !{cmd}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Step 3: Disease-related Gene Identification\n",
    "This step can identify disease related genes by analyzing patient transcriptomics' data\n",
    "\n",
    "In the `data/config_sheets` folder, create another folder called `disease`. Add an Excel file for each tissue/cell type called `disease_data_inputs_<TISSUE_NAME>`, where `<TISSUE_NAME>` is the name of the tissue you are interested in. Each sheet of this file should correspond to a separate disease to analyze using differential gene analysis. The file is formatted in the same fashion as described in the [final part of Step 1](#Importing-a-Pre-Generated-Counts-Matrix). The sheet names should be in the following format: `<DISEASE_NAME>_bulk`\n",
    "- `<DISEASE_NAME>`: This is the name of the disease you are analyzing.\n",
    "\n",
    "For example, if the disease we are interested in is lupus, and the source of the data is bulk RNA sequencing, the name of the first sheet would be `lupus_bulk`. If you are using bulk RNA sequencing, there should be a gene counts matrix file located at `data/data_matrices/<tissue_name>/<disease>` called `BulkRNAseqDataMatrix<DISEASE_NAME>_<TISSUE_NAME>`\n",
    "\n",
    "## Parameters\n",
    "- `disease_names`: The diseases you are using. This should match the first section of the sheet name in the Excel file\n",
    "- `data_source`: The datasource you are using for disease analysis. This should be`\"rnaseq\"`\n",
    "- `taxon_id`: The [NCBI Taxon ID](https://www.ncbi.nlm.nih.gov/taxonomy) to use for disease analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from como.utils import stringlist_to_list\n",
    "\n",
    "disease_names = [\"arthritis\", \"lupus_a\", \"lupus_b\"]\n",
    "data_source = \"rnaseq\"\n",
    "taxon_id = \"human\"\n",
    "\n",
    "\n",
    "for context_name in stringlist_to_list(context_names):\n",
    "    disease_config_file = f\"disease_data_inputs_{context_name}.xlsx\"\n",
    "\n",
    "    # fmt: off\n",
    "    cmd = \" \".join(\n",
    "        [\n",
    "            \"python3\", \"como/disease_analysis.py\",\n",
    "            \"--context-name\", context_name,\n",
    "            \"--config-file\", disease_config_file,\n",
    "            \"--data-source\", data_source,\n",
    "            \"--taxon-id\", str(taxon_id),\n",
    "        ]\n",
    "    )\n",
    "    # fmt: on\n",
    "\n",
    "    !{cmd}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Step 4: Drug Targets & Repurposable Drug Identification\n",
    "This step performs a series of tasks:\n",
    "1. Map drug targets in metabolic models\n",
    "2. Performs knock out simulations\n",
    "3. Compares simulation results with \"disease genes\"\n",
    "4. Identifies drug targets and repurposable drugs\n",
    "\n",
    "## Execution Steps\n",
    "### Drug Database\n",
    "A processed drug-target file is included in the `data` folder, called `Repurposing_Hub_export.txt`. If you would like to include an additional drug-target file, please model your own file after the included one. Alternatively, if you would like to update to a newer version of the database, simply export from the [Drug Repurposing Hub](https://clue.io/repurposing-app). If you do this, remove all `activators`, `agonists`, and `withdrawn` drugs. Replace the `data/Repurposing_Hub_export.txt` file.\n",
    "\n",
    "### Using Automatically Created Models\n",
    "This step will use the models generated in Step 4, above. It is **highly** recommended to use refined and validated models for further analysis (i.e., before running this step of the pipeline). If you would like to use a custom model, instead of the one created by COMO, edit the `model_files` dictionary. An example is shown here:\n",
    "```python\n",
    "model_files = {\n",
    "   \"exampleTissueModel\": \"/home/jovyan/main/data/myModels/exampleTissueModel.mat\",\n",
    "   \"anotherTissueModel\": \"/home/jovyan/main/data/myModels/anotherTissueModel.json\",\n",
    "   \"thirdTissueModel\": \"/home/jovyan/main/data/myModels/thirdTissueModel.xml\"\n",
    "}\n",
    "```\n",
    "\n",
    "❗The path `/home/jovyan/main/` **<ins>MUST</ins>** stay the same. If it does not, your model **will not be found**\n",
    "\n",
    "\n",
    "## Parameters\n",
    "Other than the `model_files` parameter (if required), the only other parameter for this section is the `solver` option\n",
    "\n",
    "- `solver`: The solver you would like to use. Available options are `\"gurobi\"` or `\"glpk\"`\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import json\n",
    "\n",
    "from como.project import Config\n",
    "from como.utils import stringlist_to_list\n",
    "from loguru import logger\n",
    "\n",
    "# Knock out simulation for the analyzed tissues and diseases\n",
    "model_files = {\n",
    "    # \"context_name\": \"/path/to/model.mat\"\n",
    "    # EXAMPLE -> \"Treg\": \"/home/jovyan/main/data/results/naiveB/naiveB_SpecificModel_IMAT.mat\"\n",
    "}\n",
    "sovler = \"gurobi\"\n",
    "\n",
    "\n",
    "config = Config()\n",
    "\n",
    "drug_raw_file = \"Repurposing_Hub_export.txt\"\n",
    "for context in stringlist_to_list(context_names):\n",
    "    for recon_algorithm in recon_algorithms:\n",
    "        for disease in disease_names:\n",
    "            disease_path = config.data_dir / \"results\" / context / disease\n",
    "            out_dir = config.data_dir / \"results\" / context / disease\n",
    "            tissue_gene_folder = config.data_dir / context\n",
    "            tissue_gene_folder.mkdir(parents=True, exist_ok=True)\n",
    "\n",
    "            if not disease_path.exists():\n",
    "                logger.warning(f\"Disease path doesn't exist! Looking for {disease_path}\")\n",
    "                continue\n",
    "\n",
    "            # load the results of step 3 to dictionary \"disease_files\"\n",
    "            step3_results_file = config.data_dir / \"results\" / context / disease / \"step2_results_files.json\"\n",
    "\n",
    "            with step3_results_file.open(\"r\") as json_file:\n",
    "                disease_files = json.load(json_file)\n",
    "                down_regulated_disease_genes = disease_files[\"down_regulated\"]\n",
    "                up_regulated_disease_genes = disease_files[\"up_regulated\"]\n",
    "\n",
    "            if context in model_files:\n",
    "                tissue_specific_model_filepath = model_files[context]\n",
    "            else:\n",
    "                tissue_specific_model_filepath = (\n",
    "                    config.data_dir / \"results\" / context / f\"{context}_SpecificModel_{recon_algorithm}.mat\"\n",
    "                )\n",
    "\n",
    "            # fmt: off\n",
    "            cmd = [\n",
    "                \"python3\", \"como/knock_out_simulation.py\",\n",
    "                \"--context-model\", tissue_specific_model_filepath,\n",
    "                \"--context-name\", context,\n",
    "                \"--disease-name\", disease,\n",
    "                \"--disease-up\", up_regulated_disease_genes,\n",
    "                \"--disease-down\", down_regulated_disease_genes,\n",
    "                \"--raw-drug-file\", drug_raw_file,\n",
    "                \"--solver\", sovler,\n",
    "                # \"--test-all\"\n",
    "            ]\n",
    "            # fmt: on\n",
    "\n",
    "            if recon_algorithm == \"IMAT\":\n",
    "                cmd.extend([\"--reference-flux-file\", config.data_dir / \"results\" / context / \"IMAT_flux.csv\"])\n",
    "\n",
    "            cmd = \" \".join(cmd)\n",
    "            !{cmd}"
   ]
  }
 ],
 "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.10.15"
  },
  "toc-autonumbering": true,
  "toc-showcode": true
 },
 "nbformat": 4,
 "nbformat_minor": 4
}