Workflow.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "a62e33fb",
   "metadata": {},
   "source": [
    "Computation Pipeline to:\n",
    "- Pre-process single cell RNA-seq dataset using Scanpy\n",
    "- Apply denoising on single cell RNA-seq data using MAGIC\n",
    "- Compute Reaction Activity Score\n",
    "- Compute single-cell Flux Balance Analysis\n",
    "- Perform Flux Cluster Analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "87bb00f5",
   "metadata": {},
   "source": [
    "Load libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "28122599",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "from scanpy import AnnData\n",
    "import scanpy as sc\n",
    "import cobra as cb\n",
    "import numpy as np\n",
    "import time\n",
    "import sys\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import os\n",
    "from scipy.sparse import csr_matrix\n",
    "import sys"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "391ffa56",
   "metadata": {},
   "source": [
    "Set denoising strategy. Three possible denoisers can be used (magic,enhance, saver). Denoising can be applied on the readcounts or on the RAS matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "44346ee6",
   "metadata": {},
   "outputs": [],
   "source": [
    "denoising_strategy=\"magic_counts\" \n",
    "compute_fluxes=True\n",
    "save_flux=True\n",
    "save_ras=True\n",
    "save_counts=False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "683bbb61",
   "metadata": {},
   "outputs": [],
   "source": [
    "if denoising_strategy==\"epsilon\":\n",
    "    eps=0.01\n",
    "else:\n",
    "    eps=0"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a5862ae3",
   "metadata": {},
   "source": [
    "Set name of the dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b14df9fc",
   "metadata": {},
   "outputs": [],
   "source": [
    "#file_input=\"datasetGSE110949\"\n",
    "file_input=\"datasetE-GEOD-86618\"\n",
    "#file_input=\"datasetGSE118056\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d27c1cc4",
   "metadata": {},
   "source": [
    "## Load RNA-seq and pre-processing before RAS computation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b0706b36",
   "metadata": {},
   "source": [
    "Load single-cell dataset and stored in AnnData object"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2245bcbd",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "adata=sc.read_h5ad(\"data/\"+file_input)\n",
    "adata"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "07d437da",
   "metadata": {},
   "source": [
    "Quality check filtering (remove cell with few expressed genes and genes not expressed in 3 cells at least)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cba97923",
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pp.filter_cells(adata, min_genes=2000)\n",
    "sc.pp.filter_genes(adata, min_cells=3)\n",
    "adata"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d68af87f",
   "metadata": {},
   "source": [
    "TPM normalization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2bc5fb45",
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pp.normalize_per_cell(adata, counts_per_cell_after=1e6)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "965e1e5a",
   "metadata": {},
   "source": [
    "Store cell labels "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e77666b9",
   "metadata": {},
   "outputs": [],
   "source": [
    "cells=list(adata.obs.index)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8b4bb9d5",
   "metadata": {},
   "source": [
    "## Application of denoising"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7f0239c1",
   "metadata": {},
   "source": [
    "Before denoising"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ed5a1238",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "adata.to_df()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "311359ff",
   "metadata": {},
   "source": [
    "Applying denoising on the readcounts matrix. \n",
    "- MAGIC can be used using the Scanpy function sc.external.pp.magic\n",
    "- To use ENHANCE, you need to download the code from the ENHANCE repository https://github.com/yanailab/enhance\n",
    "- SAVER denoiser is available as R-package. We provide an R-script to run SAVER named saver_readcounts.R"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "69a9a2bd",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "if denoising_strategy==\"magic_counts\":\n",
    "    sc.external.pp.magic(adata,random_state=0)   \n",
    "adata"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4e2e44dd",
   "metadata": {},
   "source": [
    "After denoising"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "715b8028",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "adata.to_df()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5ff10bf9",
   "metadata": {},
   "source": [
    "Save read counts data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7916accd",
   "metadata": {},
   "outputs": [],
   "source": [
    "if save_counts:\n",
    "    adata.X = csr_matrix(adata.X)\n",
    "    adata.write(\"outputs/readcounts/\"+file_input)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ee99aab0",
   "metadata": {},
   "source": [
    "## Load metabolic model"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "37d67d4a",
   "metadata": {},
   "source": [
    "Read the SBML model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b783a339",
   "metadata": {},
   "outputs": [],
   "source": [
    "if file_input==\"datasetE-GEOD-86618\":\n",
    "    model=cb.io.read_sbml_model('models/ENGRO2_ensg.xml')       #gene ensg notation\n",
    "else:\n",
    "    model=cb.io.read_sbml_model('models/ENGRO2_genesymbol.xml')       #gene ensg notation\n",
    "model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4ed2be75",
   "metadata": {},
   "outputs": [],
   "source": [
    "reactions=[reaction.id for reaction in model.reactions]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "35fd61cc",
   "metadata": {},
   "source": [
    "Set objective function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4b38528f",
   "metadata": {},
   "outputs": [],
   "source": [
    "if file_input==\"datasetE-GEOD-86618\":\n",
    "    model.objective=\"ATPM\"\n",
    "else:\n",
    "    model.objective=\"Biomass\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9dd0222c",
   "metadata": {},
   "source": [
    "## FVA computation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3a67f376",
   "metadata": {},
   "outputs": [],
   "source": [
    "from cobra.flux_analysis import flux_variability_analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "35b2b581",
   "metadata": {},
   "outputs": [],
   "source": [
    "#dfFVA=flux_variability_analysis(model,fraction_of_optimum=0)\n",
    "#dfFVA=dfFVA.round(decimals=4)\n",
    "#dfFVA[\"maxABS\"]=dfFVA.abs().max(1)\n",
    "#dfFVA.to_csv(\"data/FVA.csv\")\n",
    "#dfFVA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b292af32",
   "metadata": {},
   "outputs": [],
   "source": [
    "dfFVA=pd.read_csv(\"data/FVA.csv\",index_col=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "698efe95",
   "metadata": {},
   "outputs": [],
   "source": [
    "dfFVA"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7ff41cb8",
   "metadata": {},
   "source": [
    "## RAS computation "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d35e51d0",
   "metadata": {},
   "source": [
    "Compute Reaction Activity Scores"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ab548a2a",
   "metadata": {},
   "outputs": [],
   "source": [
    "from utils import RAS_computation as rc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8140717a",
   "metadata": {},
   "outputs": [],
   "source": [
    "adata.to_df()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "deca19be",
   "metadata": {},
   "outputs": [],
   "source": [
    "ras_object=rc(adata,model)\n",
    "ras_adata=ras_object.compute()\n",
    "ras_adata"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4f05fcd6",
   "metadata": {},
   "source": [
    "The resulting ras_adata is an anndata object where the feature now are the RAS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d89360e7",
   "metadata": {},
   "outputs": [],
   "source": [
    "ras_adata.to_df()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f5d57157",
   "metadata": {},
   "source": [
    "Applying denoising on the RAS matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0da81852",
   "metadata": {},
   "outputs": [],
   "source": [
    "if denoising_strategy==\"magic_ras\":\n",
    "    sc.external.pp.magic(ras_adata,random_state=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ea519ea1",
   "metadata": {},
   "source": [
    "After denoising"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "05f90ea5",
   "metadata": {},
   "outputs": [],
   "source": [
    "ras_adata.to_df()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "00e807cb",
   "metadata": {},
   "source": [
    "Save RAS data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3c391868",
   "metadata": {},
   "outputs": [],
   "source": [
    "if save_ras:\n",
    "    ras_adata.X = csr_matrix(ras_adata.X)\n",
    "    ras_adata.write(\"outputs/\"+file_input+\"_\"+denoising_strategy)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1db88991",
   "metadata": {},
   "source": [
    "## Single cell Flux Balance analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8a8f42d8",
   "metadata": {},
   "outputs": [],
   "source": [
    "from utils import scFBA"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "851368a6",
   "metadata": {},
   "source": [
    "Compute single cell FBA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6821d57a",
   "metadata": {},
   "outputs": [],
   "source": [
    "if compute_fluxes:\n",
    "    flux_adata=scFBA(model,ras_adata,dfFVA,eps=eps,verbose=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "601ef976",
   "metadata": {},
   "source": [
    "Save results of flux data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "661ab910",
   "metadata": {},
   "outputs": [],
   "source": [
    "if save_flux and save_flux:\n",
    "    flux_adata.write(\"outputs/fluxes/\"+file_input+\"_\"+denoising_strategy)"
   ]
  }
 ],
 "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.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}