2 - benchmark-generator.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Benchmark - Generator -\n",
    "\n",
    "## Import"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import math\n",
    "import random\n",
    "import pandas as pd\n",
    "from glob import glob\n",
    "\n",
    "import bonesis\n",
    "from functools import reduce\n",
    "\n",
    "import merrin.Loader as loader\n",
    "from merrin.Parameters import FLUXOMIC_DATA, PROTEOMIC_DATA, CINETIC_DATA, DATA_TYPE_COLUMN\n",
    "\n",
    "pd.options.mode.chained_assignment = None "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_folder(folder):\n",
    "    if not os.path.exists(folder):\n",
    "        os.makedirs(folder)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "benchmark_folder = './benchmark'\n",
    "create_folder(benchmark_folder)\n",
    "\n",
    "data_type = [\n",
    "    [CINETIC_DATA, FLUXOMIC_DATA, PROTEOMIC_DATA],\n",
    "    [CINETIC_DATA, PROTEOMIC_DATA],\n",
    "    [PROTEOMIC_DATA],\n",
    "    [CINETIC_DATA, FLUXOMIC_DATA]\n",
    "]\n",
    "\n",
    "metabolic_network_file = f'data/covert/metabolic_network.xml'\n",
    "regulatory_network_file = f'data/covert/regulatory_network.sbml'\n",
    "pkn_file = f'data/covert/interactions.txt'\n",
    "\n",
    "experiment_datasets = {\n",
    "    'covert': 'data/covert/result/pipeline_nbi/FlexFlux/out_times*.csv',\n",
    "}\n",
    "\n",
    "objective_function = 'Growth'\n",
    "hide_reactions = False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# DEGRADATION RATE\n",
    "seeds = range(0, 10)\n",
    "degradations = range(0, 60, 10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mn = loader.load_sbml(metabolic_network_file)\n",
    "inputs = sorted(mn.get_inputs())\n",
    "pkn = loader.load_pkn(pkn_file, inputs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pkn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "reactions  = mn.get_reactions()\n",
    "external_reactions = set(r for r, _ in mn.get_input_reactions())\n",
    "internal_reactions = reactions.difference(external_reactions)\n",
    "inputs = mn.get_inputs()\n",
    "outputs = mn.get_outputs()\n",
    "regulators = set(pkn.nodes()).difference(reactions).difference(inputs).difference(outputs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Search Space"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "search_space = 0\n",
    "bo = bonesis.BoNesis(bonesis.InfluenceGraph(pkn, allow_skipping_nodes=False))\n",
    "for k in inputs:\n",
    "    bo.custom(f'constant(\\\"{k}\\\", -1).')\n",
    "space = bo.local_functions().as_dict('count')\n",
    "search_space = reduce(int.__mul__, space.values())\n",
    "search_space"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data generation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def shift_external_metabolites(df):\n",
    "    to_shift = list(inputs) + ['biomass']\n",
    "    df[to_shift] = df[to_shift].shift(1)\n",
    "    df.loc[0, to_shift] = df.loc[1, to_shift]\n",
    "\n",
    "def read_simulation(csv_file: str) -> pd.DataFrame:\n",
    "    df = pd.read_csv(csv_file, sep='\\t')\n",
    "    df['Time'] = [math.floor(i * 100) if i * 100 - math.floor(i * 100)\n",
    "                  < 0.5 else math.ceil(i * 100) for i in df['Time']]\n",
    "    df.set_index('Time', inplace=True)\n",
    "    shift_external_metabolites(df)\n",
    "    return df\n",
    "\n",
    "def read_simulations(sims_path):\n",
    "    simulations = {os.path.split(os.path.splitext(f)[0])[1]: read_simulation(f)\n",
    "                   for i, f in enumerate(sorted(sims_path))}\n",
    "    return simulations\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_fluxomic(df):\n",
    "    dg = df.copy()\n",
    "    return dg[list(inputs.union(outputs).union(reactions)) + ['biomass']]\n",
    "\n",
    "def get_cinetic(df):\n",
    "    dg = df.copy()\n",
    "    return dg[list(inputs.union(outputs).union(external_reactions)) + [objective_function, 'biomass']]\n",
    "\n",
    "def get_proteomic(df):\n",
    "    threshold = 10**-5\n",
    "    dg = (df > threshold).astype(int)\n",
    "    dg = dg.drop([objective_function, 'biomass', 'Biomass'], axis=1)\n",
    "    return dg\n",
    "\n",
    "def convert_data(simulations: dict):\n",
    "    convert_inputs = {}\n",
    "    for k, sim in simulations.items():\n",
    "        fluxomic  = get_fluxomic(sim)\n",
    "        cinetic   = get_cinetic(sim)\n",
    "        proteomic = get_proteomic(sim)\n",
    "        convert_inputs[k] = (fluxomic, cinetic, proteomic)\n",
    "        \n",
    "    return convert_inputs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def extract_phases(proteomic: pd.DataFrame):\n",
    "    phases = []\n",
    "    rows = list(proteomic.iterrows())\n",
    "    current_row = (None, None)\n",
    "    for i in range(len(rows)):\n",
    "        t = i\n",
    "        state = tuple(rows[i][1].to_list())\n",
    "        if current_row[0] is None:\n",
    "            current_row = (t, state)\n",
    "        elif state != current_row[1]:\n",
    "            phases.append((current_row[0], t-1))\n",
    "            current_row = (t, state)\n",
    "    phases.append((current_row[0], len(rows)-1))\n",
    "    return phases\n",
    "\n",
    "def compress_data(simulations: dict, degradation: float, seed:int = 0):\n",
    "    random.seed(seed)\n",
    "    assert(0 <= degradation and degradation < 1)\n",
    "    compress_simulations = {}\n",
    "    for k, sim in simulations.items():\n",
    "        fluxomic, cinetic, proteomic = sim\n",
    "        phases = extract_phases(proteomic)\n",
    "\n",
    "        compress_sample = []\n",
    "        for t1, t2 in phases:\n",
    "            if t1 == t2:\n",
    "                if random.random() >= degradation:\n",
    "                    compress_sample.append(t1)\n",
    "            else:\n",
    "                n = t2 - t1 + 1\n",
    "                prob_0 = degradation**n\n",
    "                prob_1 = degradation**(n-1) * (1 - degradation)\n",
    "                prob_less_than_2 = prob_0 + prob_1\n",
    "                rand_prob = random.random()\n",
    "                if prob_less_than_2 <= rand_prob:\n",
    "                    compress_sample.extend(list(random.sample(range(t1, t2 + 1), 2)))\n",
    "                elif prob_1 <= rand_prob:\n",
    "                    compress_sample.extend(list(random.sample(range(t1, t2 + 1), 1)))\n",
    "        compress_sample.sort()\n",
    "        if compress_sample[0] != 0:\n",
    "            compress_sample = [0] + compress_sample\n",
    "        \n",
    "        c_fluxomic = fluxomic.iloc[compress_sample]\n",
    "        c_cinetic = cinetic.iloc[compress_sample]\n",
    "        c_proteomic = proteomic.iloc[compress_sample]\n",
    "\n",
    "        compress_simulations[k] = (c_fluxomic, c_cinetic, c_proteomic)\n",
    "\n",
    "    return compress_simulations\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fluxomic_noise(fluxomic: pd.DataFrame, degradation: float, seed:int = 0):\n",
    "    random.seed(seed)\n",
    "    n_fluxomic = fluxomic.applymap(lambda x: x * (1 + random.uniform(-degradation, degradation)))\n",
    "    return n_fluxomic\n",
    "\n",
    "def cinetic_noise(cinetic: pd.DataFrame, degradation: float, seed:int = 0):\n",
    "    random.seed(seed)\n",
    "    n_cinetic = cinetic.applymap(lambda x: x * (1 + random.uniform(-degradation, degradation)))\n",
    "    return n_cinetic\n",
    "\n",
    "def proteomic_noise(proteomic: pd.DataFrame, degradation: float, seed:int = 0):\n",
    "    random.seed(seed)\n",
    "    n_proteomic = proteomic.copy()\n",
    "    reactions_without_obj = reactions.union(regulators)\n",
    "    reactions_without_obj.remove(objective_function)\n",
    "    reactions_without_obj = list(reactions_without_obj)\n",
    "    n_proteomic[reactions_without_obj] = n_proteomic[reactions_without_obj].applymap(lambda x: None if random.random() < degradation else x)\n",
    "    return n_proteomic\n",
    "\n",
    "def add_noise(simulations: dict, degradation:float, seed:int = 0):\n",
    "    assert(0 <= degradation and degradation < 1)\n",
    "    random.seed(seed)\n",
    "\n",
    "    f_seed = int(random.random() * 100)\n",
    "    c_seed = int(random.random() * 100)\n",
    "    p_seed = int(random.random() * 100)\n",
    "\n",
    "    noisy_simulations = {}\n",
    "    for k, sim in simulations.items():\n",
    "        fluxomic, cinetic, proteomic = sim\n",
    "        fluxomic = fluxomic_noise(fluxomic, degradation, f_seed)\n",
    "        cinetic = cinetic_noise(cinetic, degradation, c_seed)\n",
    "        proteomic = proteomic_noise(proteomic, degradation, p_seed)\n",
    "\n",
    "        noisy_simulations[k] = (fluxomic, cinetic, proteomic)\n",
    "\n",
    "    return noisy_simulations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def drop_observations(simulations: dict, degradation:float, seed:int = 0):\n",
    "    assert(0 <= degradation and degradation < 1)\n",
    "    random.seed(seed)\n",
    "\n",
    "    drop_simulations = {}\n",
    "    for k, sim in simulations.items():\n",
    "        fluxomic, cinetic, proteomic = sim\n",
    "\n",
    "        sample = [0, 301] + random.sample(sorted(fluxomic.index[1:-1]), int((len(fluxomic.index) - 2) * (1 - degradation)))\n",
    "        sample.sort()\n",
    "\n",
    "        d_fluxomic = fluxomic.loc[sample]\n",
    "        d_cinetic = cinetic.loc[sample]\n",
    "        d_proteomic = proteomic.loc[sample]\n",
    "\n",
    "        drop_simulations[k] = (d_fluxomic, d_cinetic, d_proteomic)\n",
    "\n",
    "    return drop_simulations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def build_instance(simulations, isfluxomic:bool, iscinetic:bool, isproteomic:bool):\n",
    "    instances = {}\n",
    "    for k, sim in simulations.items():\n",
    "        fluxomic, cinetic, proteomic = sim\n",
    "\n",
    "        cinetic_col = list(inputs.union(outputs).union(external_reactions)) + [objective_function, 'biomass']\n",
    "        fluxomic_col = list(inputs.union(outputs).union(reactions)) + ['biomass']\n",
    "\n",
    "        instance = proteomic.copy()\n",
    "        if not isproteomic:\n",
    "            instance = instance.drop(list(regulators), axis=1)\n",
    "        if isfluxomic:\n",
    "            instance[fluxomic_col] = fluxomic[fluxomic_col]\n",
    "        if iscinetic or isfluxomic:\n",
    "            instance[cinetic_col] = cinetic[cinetic_col]\n",
    "\n",
    "\n",
    "        tags = []\n",
    "        if isfluxomic: tags.append(FLUXOMIC_DATA)\n",
    "        if iscinetic: tags.append(CINETIC_DATA)\n",
    "        if isproteomic: tags.append(PROTEOMIC_DATA)\n",
    "        instance[DATA_TYPE_COLUMN] = ','.join(tags)\n",
    "\n",
    "        instances[k] = instance\n",
    "\n",
    "    return instances"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Benchmark generation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ref_simulations = read_simulations(glob(experiment_datasets['covert']))\n",
    "ref_convert_simulations = convert_data(ref_simulations)\n",
    "\n",
    "for seed in seeds:\n",
    "    for degradation in degradations:\n",
    "        degradation /= 100\n",
    "        compress_simulations = compress_data(ref_convert_simulations, degradation, seed)\n",
    "        noise_simulations = add_noise(compress_simulations, degradation, seed)\n",
    "        for tags in data_type:\n",
    "            instance = build_instance(\n",
    "                noise_simulations,\n",
    "                FLUXOMIC_DATA in tags,\n",
    "                CINETIC_DATA in tags,\n",
    "                PROTEOMIC_DATA in tags\n",
    "            )\n",
    "            tags.sort()\n",
    "            tags_str = ''.join([t[0] for t in tags])\n",
    "            export_path = benchmark_folder + f'/instances/datatype-{tags_str}/degradation-{int(degradation*100)}%/seed-{seed}/'\n",
    "            create_folder(export_path)\n",
    "\n",
    "            for k, sim in instance.items():\n",
    "                sim.to_csv(export_path + k + '.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "d4d1e4263499bec80672ea0156c357c1ee493ec2b1c70f0acce89fc37c4a6abe"
  },
  "kernelspec": {
   "display_name": "Python 3",
   "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.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}