Composing models/Merge logical models.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "dc95313b-b406-4b8d-8fbb-1d2bfe56f335",
   "metadata": {},
   "source": [
    "# Merging Boolean Networks\n",
    "\n",
    "## Overview\n",
    "\n",
    "This Jupyter notebook provides tools for merging multiple Boolean network models using  `OR`, `AND`, and `Inhibitor Wins`. The notebook includes functions for:\n",
    "\n",
    "- **Parsing Boolean Expressions:** Extracting activators and inhibitors from Boolean expressions.\n",
    "- **Merging Boolean Networks:** Combining Boolean network models using different methods (`OR`, `AND`, and `Inhibitor Wins`).\n",
    "- **Writing Merged Networks to File:** Saving the final merged network to a text file or SBML-qual file.\n",
    "\n",
    "## Usage\n",
    "\n",
    "1. **Input format:**\n",
    "   - Text file similar to the Python package `Boolnet`, as described [here](https://rdrr.io/cran/BoolNet/man/loadNetwork.html)\n",
    "   - SBML-qual format\n",
    "1. **Define Input Networks:** Provide the Boolean networks as dictionaries where keys are gene names and values are Boolean expressions.\n",
    "2. **Merge Networks:** Use the `merge_networks` function to merge the networks using the desired method.\n",
    "3. **Save to File:** Use the `write_network_to_file` function to save the merged network.\n",
    "   - Text file similar to input\n",
    "   - SBML-qual format\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "0c312f06-e281-44c7-aba3-37746422a694",
   "metadata": {},
   "outputs": [],
   "source": [
    "home = \"persistent/\" # home path for Logic Model Merger\n",
    "model1name = \"Palma2021\"\n",
    "model2name = \"Ikonomi2020\"\n",
    "sbml_file1 = home + \"LogicModelMerger/Models/\" + model1name + \".sbml\"\n",
    "sbml_file2 = home + \"LogicModelMerger/Models/\" + model2name + \".sbml\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "5c940331",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Need to install libsbml in Colomoto environment\n",
    "#!conda install -c conda-forge python-libsbml"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "ba1ab538-4550-43a7-8ca7-2b05cc4dd73d",
   "metadata": {},
   "outputs": [],
   "source": [
    "import re\n",
    "import libsbml\n",
    "import os\n",
    "import boolean\n",
    "from rpy2.robjects.packages import importr\n",
    "boolnet = importr(\"BoolNet\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "4114f1b1-4bb1-46cb-b88a-02a9e8f2d907",
   "metadata": {},
   "outputs": [],
   "source": [
    "def read_network(file_path, body_separator=\",\"):\n",
    "    \"\"\"\n",
    "    Reads a Boolean network model from a text file or an SBML-qual file and returns a dictionary representation.\n",
    "    The Boolean expressions are cleaned of unnecessary parentheses.\n",
    "    \n",
    "    :param file_path: Path to the model file (text file or SBML-qual file).\n",
    "    :param body_separator: Character used to separate targets and factors in the text file.\n",
    "    :return: Dictionary with gene names as keys and their cleaned Boolean expressions as values.\n",
    "    \"\"\"\n",
    "    # Determine the file type based on the extension\n",
    "    file_extension = os.path.splitext(file_path)[1].lower()\n",
    "    \n",
    "    if file_extension in [\".txt\", \".csv\"]:\n",
    "        return read_text_network(file_path, body_separator)\n",
    "    elif file_extension in [\".xml\", \".sbml\"]:\n",
    "        return read_sbml_qual_network(file_path)\n",
    "    else:\n",
    "        raise ValueError(\"Unsupported file format. Please provide a .txt, .csv, .xml, or .sbml file.\")\n",
    "\n",
    "def read_text_network(file_path, body_separator=\",\"):\n",
    "    \"\"\"\n",
    "    Reads a Boolean network model from a text file and returns a dictionary representation.\n",
    "    The Boolean expressions are cleaned of unnecessary parentheses.\n",
    "    \n",
    "    :param file_path: Path to the model file.\n",
    "    :param body_separator: Character used to separate targets and factors in the file.\n",
    "    :return: Dictionary with gene names as keys and their cleaned Boolean expressions as values.\n",
    "    \"\"\"\n",
    "    network = {}\n",
    "\n",
    "    with open(file_path, 'r') as file:\n",
    "        lines = file.readlines()\n",
    "        \n",
    "    # Strip comments and empty lines\n",
    "    lines = [re.sub(r'#.*', '', line).strip() for line in lines if line.strip() and not line.startswith(\"#\")]\n",
    "    \n",
    "    # Check if the file is empty after stripping comments and empty lines\n",
    "    if not lines:\n",
    "        raise ValueError(\"The file is empty or contains only comments.\")\n",
    "    \n",
    "    # Check and parse the header\n",
    "    header = [part.strip() for part in lines[0].split(body_separator)]\n",
    "    if len(header) < 2 or header[0].lower() != \"targets\" or header[1].lower() not in [\"functions\", \"factors\"]:\n",
    "        raise ValueError(f\"Invalid header: {lines[0]}\")\n",
    "    \n",
    "    # Process each rule\n",
    "    for line in lines[1:]:\n",
    "        # Split the line based on the body_separator, but only at the first occurrence\n",
    "        parts = [part.strip() for part in line.split(body_separator, 1)]\n",
    "        if len(parts) < 2:\n",
    "            raise ValueError(f\"Invalid rule format: {line}\")\n",
    "        \n",
    "        target = parts[0].strip()\n",
    "        expression = parts[1].strip()\n",
    "        \n",
    "        # Validate gene names\n",
    "        if not re.match(r'^[a-zA-Z_][a-zA-Z0-9_]*$', target):\n",
    "            raise ValueError(f\"Invalid gene name: {target}\")\n",
    "        \n",
    "        network[target] = simplify_expression(expression)\n",
    "    \n",
    "    # Capitalize all the genes and expressions to make the networks case-insensitive\n",
    "    network = {gene.upper(): expression.upper() for gene, expression in network.items()}\n",
    "    \n",
    "    return network\n",
    "\n",
    "def read_sbml_qual_network(file_path):\n",
    "    \"\"\"\n",
    "    Reads a Boolean network model from an SBML-qual file and returns a dictionary representation.\n",
    "    The Boolean expressions are cleaned of unnecessary parentheses.\n",
    "    \n",
    "    :param file_path: Path to the SBML-qual file.\n",
    "    :return: Dictionary with gene names as keys and their cleaned Boolean expressions as values.\n",
    "    \"\"\"\n",
    "    reader = libsbml.SBMLReader()\n",
    "    document = reader.readSBML(file_path)\n",
    "    model = document.getModel()\n",
    "\n",
    "    # Parse qualitative species (genes)\n",
    "    genes = {}\n",
    "    listOfSpecies = model.getPlugin(\"qual\").getListOfQualitativeSpecies()\n",
    "\n",
    "    for i in range(listOfSpecies.size()):\n",
    "        species = listOfSpecies.get(i)\n",
    "        gene_name = species.getId()\n",
    "        genes[gene_name] = species.getName() or gene_name\n",
    "\n",
    "    # Parse transitions (interactions)\n",
    "    network = {}\n",
    "    listOfTransitions = model.getPlugin(\"qual\").getListOfTransitions()\n",
    "\n",
    "    for i in range(listOfTransitions.size()):\n",
    "        transition = listOfTransitions.get(i)\n",
    "        output_list = transition.getListOfOutputs()\n",
    "        function_terms = transition.getListOfFunctionTerms()\n",
    "\n",
    "        # Get function terms and convert to Boolean expressions\n",
    "        expression = []\n",
    "        for j in range(function_terms.size()):\n",
    "            function_term = function_terms.get(j)\n",
    "            math = function_term.getMath()\n",
    "            expression.append(parse_mathml(math, genes))\n",
    "        # print(expression)\n",
    "        if len(expression) > 1:\n",
    "            expression = \" | \".join(expression)\n",
    "        else:\n",
    "            expression = expression[0]\n",
    "\n",
    "        for output in output_list:\n",
    "            gene = genes[output.getQualitativeSpecies()]\n",
    "            network[gene] = simplify_expression(expression)\n",
    "\n",
    "    # Capitalize all the genes and expressions to make the networks case-insensitive\n",
    "    network = {gene.upper(): expression.upper() for gene, expression in network.items()}\n",
    "    return network\n",
    "\n",
    "\n",
    "def parse_mathml(math, genes):\n",
    "    \"\"\"\n",
    "    Parses MathML expressions to Boolean expressions.\n",
    "    \n",
    "    :param math: MathML object from libSBML.\n",
    "    :param genes: Dictionary of gene names.\n",
    "    :return: Boolean expression as a string.\n",
    "    \"\"\"\n",
    "    name = math.getType()\n",
    "    \n",
    "    if name == libsbml.AST_CONSTANT_TRUE:\n",
    "        return \"1\"\n",
    "    elif name == libsbml.AST_CONSTANT_FALSE:\n",
    "        return \"0\"\n",
    "    elif name == libsbml.AST_NAME:\n",
    "        return genes.get(math.getName(), math.getName())\n",
    "    elif name == libsbml.AST_REAL:\n",
    "        return str(math.getReal())\n",
    "    elif name == libsbml.AST_INTEGER:\n",
    "        return str(math.getInteger())\n",
    "    elif name == libsbml.AST_RELATIONAL_EQ:\n",
    "        left_child = parse_mathml(math.getChild(0), genes)\n",
    "        right_child = parse_mathml(math.getChild(1), genes)\n",
    "        if right_child == \"1\":\n",
    "            return left_child  # GATA1 == 1 should return GATA1\n",
    "        elif right_child == \"0\":\n",
    "            return f\"!{left_child}\"  # GATA1 == 0 should return !GATA1\n",
    "        else:\n",
    "            return f\"({left_child} == {right_child})\"\n",
    "    elif name == libsbml.AST_LOGICAL_AND:\n",
    "        children = [parse_mathml(math.getChild(i), genes) for i in range(math.getNumChildren())]\n",
    "        return f\"({' & '.join(children)})\"\n",
    "    elif name == libsbml.AST_LOGICAL_OR:\n",
    "        children = [parse_mathml(math.getChild(i), genes) for i in range(math.getNumChildren())]\n",
    "        return f\"({' | '.join(children)})\"\n",
    "    elif name == libsbml.AST_LOGICAL_NOT:\n",
    "        child = parse_mathml(math.getChild(0), genes)\n",
    "        if math.getChild(0).getType() in (libsbml.AST_LOGICAL_AND, libsbml.AST_LOGICAL_OR):\n",
    "            return f\"!({child})\"\n",
    "        else:\n",
    "            return f\"!{child}\"\n",
    "    else:\n",
    "        raise ValueError(f\"Unsupported MathML element: {libsbml.formulaToString(math)}\")\n",
    "\n",
    "\n",
    "def simplify_expression(expression):\n",
    "    \"\"\"\n",
    "    Simplifies the Boolean expression by removing redundant parentheses using boolean.py package.\n",
    "    \n",
    "    :param expression: The Boolean expression as a string.\n",
    "    :return: The simplified expression.\n",
    "    \"\"\"\n",
    "    algebra = boolean.BooleanAlgebra()\n",
    "    algebra.parse(expression)\n",
    "    simplified_expression = str(algebra.parse(expression).simplify())\n",
    "    # Replace '~' with '!'\n",
    "    simplified_expression = simplified_expression.replace('~', '!')\n",
    "\n",
    "    return simplified_expression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "5301dcfc-321e-44d4-a800-2f9d3d1f7803",
   "metadata": {},
   "outputs": [],
   "source": [
    "def parse_expression(expression):\n",
    "    \"\"\"\n",
    "    Parses a Boolean expression to separate activators and inhibitors.\n",
    "    \"\"\"\n",
    "    activators = []\n",
    "    inhibitors = []\n",
    "\n",
    "    # Detect complex inhibitors like !(GATA1 & GATA2)\n",
    "    complex_inhibitors = re.findall(r'!\\(([^)]+)\\)', expression)\n",
    "    for complex_inhibitor in complex_inhibitors:\n",
    "        inhibitors.append(complex_inhibitor.strip())\n",
    "        expression = expression.replace(f\"!({complex_inhibitor.strip()})\", '')\n",
    "\n",
    "    # Detect simple inhibitors like !GATA1\n",
    "    simple_inhibitors = re.findall(r'![a-zA-Z_][a-zA-Z0-9_]*', expression)\n",
    "    inhibitors.extend([inhibitor[1:] for inhibitor in simple_inhibitors])\n",
    "\n",
    "    # Remove all inhibitors from the expression to identify activators\n",
    "    for inhibitor in inhibitors:\n",
    "        expression = expression.replace(f\"&!{inhibitor}\", '').replace(f\"|!{inhibitor}\", '').replace(f\"!{inhibitor}\", '')\n",
    "\n",
    "    # Detect activator groups connected by '&' or within parentheses\n",
    "    activator_groups = re.findall(r'\\([^()]*\\)|[^()]*&[^()]*', expression)\n",
    "    # Remove the symbol & behind the last activator name and filter out empty strings and exactly '(|)'\n",
    "    activator_groups = [group.rstrip('&').strip() for group in activator_groups if group.strip() != '(|)']\n",
    "    \n",
    "    for group in activator_groups:\n",
    "        if group.strip():\n",
    "            activators.append(group.strip())\n",
    "        expression = expression.replace(group, '')\n",
    "\n",
    "    # Detect remaining simple activators\n",
    "    activators.extend(re.findall(r'[a-zA-Z_][a-zA-Z0-9_]*', expression))\n",
    "\n",
    "    # Remove empty parentheses from the activators list\n",
    "    activators = [activator for activator in activators if activator != '()']\n",
    "    # Remove symbol if it is the first character\n",
    "    activators = [activator[1:] if activator[0] in ['&', '|'] else activator for activator in activators]\n",
    "    # Remove symbol and the parentheses if it is the first character within parentheses\n",
    "    activators = [activator[2:-1] if activator.startswith('(&') or activator.startswith('(|') else activator for activator in activators]\n",
    "    # Remove duplicates\n",
    "    activators = list(set(activators))\n",
    "    inhibitors = list(set(inhibitors))\n",
    "    return activators, inhibitors\n",
    "\n",
    "def check_inhibitor_wins_rule(expression):\n",
    "    \"\"\"\n",
    "    Ensures that the \"Inhibitor Wins\" rule is respected within a single expression.\n",
    "    \"\"\"\n",
    "    components = re.split(r'\\|', expression)\n",
    "    \n",
    "    for component in components:\n",
    "        component = component.strip()\n",
    "        if '&' in component:\n",
    "            activators, inhibitors = parse_expression(component)\n",
    "            if activators and inhibitors:\n",
    "                continue\n",
    "        if \"!\" not in component:\n",
    "            continue\n",
    "        activators, inhibitors = parse_expression(component)\n",
    "        if activators and inhibitors:\n",
    "            return False\n",
    "    \n",
    "    return True\n",
    "\n",
    "def merge_networks(networks, method=\"OR\", descriptive=True):\n",
    "    \"\"\"\n",
    "    Merges multiple Boolean network models using the specified method (OR, AND, or Inhibitor Wins).\n",
    "    Outputs a description of the merged model if descriptive is True.\n",
    "    \"\"\"\n",
    "    if method not in [\"OR\", \"AND\", \"Inhibitor Wins\"]:\n",
    "        raise ValueError(\"Invalid method. Use 'OR', 'AND', or 'Inhibitor Wins'.\")\n",
    "\n",
    "    merged_network = {}\n",
    "    overlap_genes = set()\n",
    "    descriptions = {}\n",
    "    warning_gene = {}\n",
    "    algebra = boolean.BooleanAlgebra()\n",
    "\n",
    "    # Track the number of genes in each individual model\n",
    "    individual_gene_counts = [len(network) for network in networks]\n",
    "\n",
    "    # Iterate through each network\n",
    "    for idx, network in enumerate(networks):\n",
    "        for gene, expression in network.items():\n",
    "            if gene in merged_network:\n",
    "                overlap_genes.add(gene)\n",
    "                # Ignore if the expression is the same\n",
    "                if algebra.parse(merged_network[gene]) == algebra.parse(expression):\n",
    "                    if descriptive:\n",
    "                        descriptions.setdefault(gene, []).append((f\"Model {idx + 1}\", expression))\n",
    "                        descriptions[gene].append((\"Merged\", merged_network[gene]))\n",
    "                    continue\n",
    "                # Determine how to merge the expressions based on the method\n",
    "                if method == \"OR\":\n",
    "                    merged_expression = f\"({merged_network[gene]})|({expression})\"\n",
    "                elif method == \"AND\":\n",
    "                    merged_expression = f\"({merged_network[gene]})&({expression})\"\n",
    "                elif method == \"Inhibitor Wins\":\n",
    "                    activators, inhibitors = parse_expression(expression)\n",
    "                    existing_activators, existing_inhibitors = parse_expression(merged_network[gene])\n",
    "                    combined_activators = list(set(existing_activators + activators))\n",
    "                    combined_inhibitors = list(set(existing_inhibitors + inhibitors))\n",
    "                    \n",
    "                    activator_expr = '|'.join(filter(None, combined_activators)) if combined_activators else \"\"\n",
    "                    inhibitor_expr = '|'.join(filter(None, combined_inhibitors)) if combined_inhibitors else \"\"\n",
    "                    \n",
    "                    activator_genes = set(re.split(r'\\|', activator_expr))\n",
    "                    inhibitor_genes = set(re.split(r'\\|', inhibitor_expr))\n",
    "                    if activator_genes & inhibitor_genes:\n",
    "                        overlappings = list(activator_genes & inhibitor_genes)\n",
    "                        for overlapping in overlappings:\n",
    "                            activator_expr = re.sub(rf'(\\||&)?{overlapping}(\\||&)?', '', activator_expr).strip('|&')\n",
    "                        warning_gene[gene] = overlappings\n",
    "                    \n",
    "                    if inhibitor_expr:\n",
    "                        if '&' in inhibitor_expr or '|' in inhibitor_expr:\n",
    "                            merged_expression = f\"({activator_expr})&!({inhibitor_expr})\"\n",
    "                        else:\n",
    "                            merged_expression = f\"({activator_expr})&!{inhibitor_expr}\"\n",
    "                    else:\n",
    "                        merged_expression = activator_expr\n",
    "                else:\n",
    "                    merged_expression = expression\n",
    "                \n",
    "                # Simplify the merged expression\n",
    "                # print(gene)\n",
    "                # print(merged_expression)\n",
    "                merged_network[gene] = simplify_expression(merged_expression)\n",
    "\n",
    "                if descriptive:\n",
    "                    descriptions.setdefault(gene, []).append((f\"Model {idx + 1}\", expression))\n",
    "                    descriptions[gene].append((\"Merged\", merged_network[gene]))\n",
    "                \n",
    "            else:\n",
    "                merged_network[gene] = expression\n",
    "                if descriptive:\n",
    "                    descriptions[gene] = [(f\"Model {idx + 1}\", expression)]\n",
    "\n",
    "    if descriptive:\n",
    "        print(f\"Merging Method: {method}\")\n",
    "        print(f\"Total Genes in Merged Network: {len(merged_network)}\")\n",
    "        print(f\"Number of Genes in Each Individual Model:\")\n",
    "        for i, count in enumerate(individual_gene_counts, 1):\n",
    "            print(f\"  Model {i}: {count} genes\")\n",
    "        print(f\"Overlapping Genes: {len(overlap_genes)}\")\n",
    "        if overlap_genes:\n",
    "            print(f\"Overlapping Genes List: {', '.join(overlap_genes)}\")\n",
    "            for gene in overlap_genes:\n",
    "                print(f\"\\nGene: {gene}\")\n",
    "                if gene in warning_gene.keys():\n",
    "                    print(f\"Warning: possible conflicts for {warning_gene[gene]}, keeping only inhibitor.\")\n",
    "                for desc in descriptions[gene]:\n",
    "                    print(f\"  {desc[0]} Function: {desc[1]}\")\n",
    "        else:\n",
    "            print(\"No overlapping genes found.\")\n",
    "     \n",
    "    return merged_network"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "49dcf21a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Network #1:\n",
      "{'FLT3': 'FLT3', 'NPM1': 'NPM1', 'AKT1': 'FLT3', 'APOPTOSIS': '!BCL2&TP53', 'BCL2': 'MAPK1&!TP53', 'CCND1': '!(DNMT3A|GSK3B)', 'CDKN2A': 'NPM1', 'CEBPA': '!FLT3', 'DIFFERENTIATION': '!MEIS1&(CEBPA|ETV6)', 'DNMT3A': 'DNMT3A', 'MAPK1': 'FLT3', 'ETV6': '!MAPK1', 'FBXW7': 'NPM1', 'GSK3B': '!AKT1', 'HOXA9': '!NPM1', 'MEIS1': '!DNMT3A&HOXA9', 'MYC': 'MAPK1&(!FBXW7|!GSK3B)', 'PROLIFERATION': 'CCND1|MEIS1|MYC|SOX4|STAT5A', 'SOX4': '!CEBPA', 'STAT5A': 'FLT3', 'TP53': 'CDKN2A'}\n",
      "Network #2:\n",
      "{'EXTERNAL_QUIESCENCE': 'EXTERNAL_QUIESCENCE', 'EXTERNAL_CYCLING': 'EXTERNAL_CYCLING', 'PIK3CA': 'NRAS', 'TSC1': '!AKT1', 'MTOR': '!TSC1', 'FOXO3': '!AKT1&EXTERNAL_QUIESCENCE', 'ATM': 'FOXO3', 'ROS1': '!ATM|!BMI1|!FOXO3|MITOCHONDRIA|!TP53', 'MITOCHONDRIA': 'MTOR', 'AUTOPHAGY': 'FOXO3&!MTOR&ROS1', 'NRAS': 'EXTERNAL_CYCLING', 'ETS1': '!MEFV&NRAS', 'MEFV': 'NRAS', 'GSK3B': '!AKT1', 'CTNNB1': '!GSK3B', 'MYC': 'CTNNB1&!GSK3B', 'BMI1': 'MYC|(ATM&FOXO3)', 'MDM2': '!ATM&!CDKN2D&(MEFV|TP53)', 'TP53': '!MDM2', 'CDKN1C': 'EXTERNAL_QUIESCENCE|FOXO3', 'CDKN1A': '!MYC&(EXTERNAL_QUIESCENCE|FOXO3|GFI1|TP53)', 'CDKN1B': 'FOXO3', 'GFI1': 'TP53', 'RB1': '!CCND1&!CCNE1', 'E2F1': '!GFI1&!RB1', 'CCND1': '!CDKN1C&!CDKN2A&MYC', 'CCNE1': '!CDKN1C&E2F1&(CCND1|(!CDKN1A&!CDKN1B))', 'S_PHASE': 'CCNE1&E2F1', 'AKT1': 'PIK3CA', 'CDKN2D': '!BMI1&(E2F1|ROS1)', 'CDKN2A': '!BMI1&(ETS1|ROS1)', 'PRO_APOPTOTIC_PROTEINS': '!AKT1&ROS1&TP53', 'ANTI_APOPTOTIC_PROTEINS': '!GSK3B&(EXTERNAL_QUIESCENCE|NRAS)', 'CYCS': '!ANTI_APOPTOTIC_PROTEINS&PRO_APOPTOTIC_PROTEINS', 'APOPTOSIS': '!AKT1&CYCS', 'SENESCENCE': '(CDKN1A&ROS1&TP53)|(CDKN2A&ROS1)'}\n"
     ]
    }
   ],
   "source": [
    "# Load the networks\n",
    "network1 = read_network(sbml_file1)\n",
    "network2 = read_network(sbml_file2)\n",
    "print(\"Network #1:\")\n",
    "print(network1)\n",
    "print(\"Network #2:\")\n",
    "print(network2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "fc8b75b2-d5ef-41ff-a177-577659646c99",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Merging Method: OR\n",
      "Total Genes in Merged Network: 50\n",
      "Number of Genes in Each Individual Model:\n",
      "  Model 1: 21 genes\n",
      "  Model 2: 36 genes\n",
      "Overlapping Genes: 7\n",
      "Overlapping Genes List: GSK3B, TP53, MYC, CDKN2A, APOPTOSIS, AKT1, CCND1\n",
      "\n",
      "Gene: GSK3B\n",
      "  Model 1 Function: !AKT1\n",
      "  Model 2 Function: !AKT1\n",
      "  Merged Function: !AKT1\n",
      "\n",
      "Gene: TP53\n",
      "  Model 1 Function: CDKN2A\n",
      "  Model 2 Function: !MDM2\n",
      "  Merged Function: CDKN2A|!MDM2\n",
      "\n",
      "Gene: MYC\n",
      "  Model 1 Function: MAPK1&(!FBXW7|!GSK3B)\n",
      "  Model 2 Function: CTNNB1&!GSK3B\n",
      "  Merged Function: (CTNNB1&!GSK3B)|(MAPK1&(!FBXW7|!GSK3B))\n",
      "\n",
      "Gene: CDKN2A\n",
      "  Model 1 Function: NPM1\n",
      "  Model 2 Function: !BMI1&(ETS1|ROS1)\n",
      "  Merged Function: NPM1|(!BMI1&(ETS1|ROS1))\n",
      "\n",
      "Gene: APOPTOSIS\n",
      "  Model 1 Function: !BCL2&TP53\n",
      "  Model 2 Function: !AKT1&CYCS\n",
      "  Merged Function: (!AKT1&CYCS)|(!BCL2&TP53)\n",
      "\n",
      "Gene: AKT1\n",
      "  Model 1 Function: FLT3\n",
      "  Model 2 Function: PIK3CA\n",
      "  Merged Function: FLT3|PIK3CA\n",
      "\n",
      "Gene: CCND1\n",
      "  Model 1 Function: !(DNMT3A|GSK3B)\n",
      "  Model 2 Function: !CDKN1C&!CDKN2A&MYC\n",
      "  Merged Function: (!CDKN1C&!CDKN2A&MYC)|(!DNMT3A&!GSK3B)\n"
     ]
    }
   ],
   "source": [
    "# Example usage:\n",
    "\n",
    "# Merge using OR approach\n",
    "merged_network_or = merge_networks([network1, network2], method=\"OR\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "28d97a66-8f73-4eab-ada0-e9094f1b1b31",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Merging Method: AND\n",
      "Total Genes in Merged Network: 50\n",
      "Number of Genes in Each Individual Model:\n",
      "  Model 1: 21 genes\n",
      "  Model 2: 36 genes\n",
      "Overlapping Genes: 7\n",
      "Overlapping Genes List: GSK3B, TP53, MYC, CDKN2A, APOPTOSIS, AKT1, CCND1\n",
      "\n",
      "Gene: GSK3B\n",
      "  Model 1 Function: !AKT1\n",
      "  Model 2 Function: !AKT1\n",
      "  Merged Function: !AKT1\n",
      "\n",
      "Gene: TP53\n",
      "  Model 1 Function: CDKN2A\n",
      "  Model 2 Function: !MDM2\n",
      "  Merged Function: CDKN2A&!MDM2\n",
      "\n",
      "Gene: MYC\n",
      "  Model 1 Function: MAPK1&(!FBXW7|!GSK3B)\n",
      "  Model 2 Function: CTNNB1&!GSK3B\n",
      "  Merged Function: CTNNB1&!GSK3B&MAPK1\n",
      "\n",
      "Gene: CDKN2A\n",
      "  Model 1 Function: NPM1\n",
      "  Model 2 Function: !BMI1&(ETS1|ROS1)\n",
      "  Merged Function: !BMI1&NPM1&(ETS1|ROS1)\n",
      "\n",
      "Gene: APOPTOSIS\n",
      "  Model 1 Function: !BCL2&TP53\n",
      "  Model 2 Function: !AKT1&CYCS\n",
      "  Merged Function: !AKT1&!BCL2&CYCS&TP53\n",
      "\n",
      "Gene: AKT1\n",
      "  Model 1 Function: FLT3\n",
      "  Model 2 Function: PIK3CA\n",
      "  Merged Function: FLT3&PIK3CA\n",
      "\n",
      "Gene: CCND1\n",
      "  Model 1 Function: !(DNMT3A|GSK3B)\n",
      "  Model 2 Function: !CDKN1C&!CDKN2A&MYC\n",
      "  Merged Function: !CDKN1C&!CDKN2A&!DNMT3A&!GSK3B&MYC\n"
     ]
    }
   ],
   "source": [
    "# Merge using AND approach\n",
    "merged_network_and = merge_networks([network1, network2], method=\"AND\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "e4680cf1-1072-4231-8d49-9a55d002b27d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Merging Method: Inhibitor Wins\n",
      "Total Genes in Merged Network: 50\n",
      "Number of Genes in Each Individual Model:\n",
      "  Model 1: 21 genes\n",
      "  Model 2: 36 genes\n",
      "Overlapping Genes: 7\n",
      "Overlapping Genes List: GSK3B, TP53, MYC, CDKN2A, APOPTOSIS, AKT1, CCND1\n",
      "\n",
      "Gene: GSK3B\n",
      "  Model 1 Function: !AKT1\n",
      "  Model 2 Function: !AKT1\n",
      "  Merged Function: !AKT1\n",
      "\n",
      "Gene: TP53\n",
      "  Model 1 Function: CDKN2A\n",
      "  Model 2 Function: !MDM2\n",
      "  Merged Function: CDKN2A&!MDM2\n",
      "\n",
      "Gene: MYC\n",
      "  Model 1 Function: MAPK1&(!FBXW7|!GSK3B)\n",
      "  Model 2 Function: CTNNB1&!GSK3B\n",
      "  Merged Function: !FBXW7&!GSK3B&(CTNNB1|MAPK1)\n",
      "\n",
      "Gene: CDKN2A\n",
      "  Model 1 Function: NPM1\n",
      "  Model 2 Function: !BMI1&(ETS1|ROS1)\n",
      "  Merged Function: !BMI1&(ETS1|NPM1|ROS1)\n",
      "\n",
      "Gene: APOPTOSIS\n",
      "  Model 1 Function: !BCL2&TP53\n",
      "  Model 2 Function: !AKT1&CYCS\n",
      "  Merged Function: !AKT1&!BCL2&(CYCS|TP53)\n",
      "\n",
      "Gene: AKT1\n",
      "  Model 1 Function: FLT3\n",
      "  Model 2 Function: PIK3CA\n",
      "  Merged Function: FLT3|PIK3CA\n",
      "\n",
      "Gene: CCND1\n",
      "  Model 1 Function: !(DNMT3A|GSK3B)\n",
      "  Model 2 Function: !CDKN1C&!CDKN2A&MYC\n",
      "  Merged Function: !CDKN1C&!CDKN2A&!DNMT3A&!GSK3B&MYC\n"
     ]
    }
   ],
   "source": [
    "# Merge using Inhibitor Wins approach\n",
    "merged_network_inhibitor_wins = merge_networks([network1, network2], method=\"Inhibitor Wins\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "453520b2-ab60-44b6-bbad-171cb0822649",
   "metadata": {},
   "outputs": [],
   "source": [
    "def write_network_to_file(network, filename, format=\"text\"):\n",
    "    \"\"\"\n",
    "    Writes the merged Boolean network to a file in the specified format.\n",
    "    \n",
    "    :param network: Dictionary representing the merged Boolean network.\n",
    "    :param filename: The name of the file to write the network to.\n",
    "    :param format: The format of the output file, either \"text\" or \"sbml\".\n",
    "    \"\"\"\n",
    "    if format == \"text\":\n",
    "        with open(filename + \".txt\", 'w') as file:\n",
    "            file.write(\"targets, factors\\n\")\n",
    "            for gene, expression in network.items():\n",
    "                file.write(f\"{gene}, {expression}\\n\")\n",
    "        print(f\"Network successfully written to {filename}.txt\")\n",
    "    \n",
    "    elif format == \"sbml\":\n",
    "        # Write to a temporary text file\n",
    "        temp_text_file = \"temp_network.txt\"\n",
    "        with open(temp_text_file, 'w') as file:\n",
    "            file.write(\"targets, factors\\n\")\n",
    "            for gene, expression in network.items():\n",
    "                file.write(f\"{gene}, {expression}\\n\")\n",
    "        \n",
    "        # Load the network from the text file using boolnet\n",
    "        net = boolnet.loadNetwork(temp_text_file)\n",
    "        \n",
    "        # Export to SBML format\n",
    "        boolnet.toSBML(net, filename + \".sbml\")\n",
    "        \n",
    "        # Delete the temporary text file\n",
    "        os.remove(temp_text_file)\n",
    "        \n",
    "        print(f\"Network successfully written to {filename}.sbml\")\n",
    "    \n",
    "    else:\n",
    "        raise ValueError(\"Invalid format specified. Use 'text' or 'sbml'.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "5a20c879-19b3-4853-bf46-c8279cb2e0f6",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Network successfully written to ../Models/merged_and_Palma2021_Ikonomi2020.txt\n",
      "Network successfully written to ../Models/merged_inhibitor_wins_Palma2021_Ikonomi2020.txt\n",
      "Network successfully written to ../Models/merged_or_Palma2021_Ikonomi2020.txt\n",
      "Network successfully written to ../Models/merged_and_Palma2021_Ikonomi2020.sbml\n",
      "Network successfully written to ../Models/merged_inhibitor_wins_Palma2021_Ikonomi2020.sbml\n",
      "Network successfully written to ../Models/merged_or_Palma2021_Ikonomi2020.sbml\n"
     ]
    }
   ],
   "source": [
    "merged_and_name = \"merged_and_\" + model1name + \"_\" + model2name\n",
    "merged_inhibitor_wins_name = \"merged_inhibitor_wins_\" + model1name + \"_\" + model2name\n",
    "merged_or_name = \"merged_or_\" + model1name + \"_\" + model2name\n",
    "\n",
    "write_network_to_file(merged_network_and, home + \"LogicModelMerger/Models/\" + merged_and_name, format=\"text\")\n",
    "write_network_to_file(merged_network_inhibitor_wins, home + \"LogicModelMerger/Models/\" + merged_inhibitor_wins_name, format=\"text\")\n",
    "write_network_to_file(merged_network_or, home + \"LogicModelMerger/Models/\" + merged_or_name, format=\"text\")\n",
    "write_network_to_file(merged_network_and, home + \"LogicModelMerger/Models/\" + merged_and_name, format=\"sbml\")\n",
    "write_network_to_file(merged_network_inhibitor_wins, home + \"LogicModelMerger/Models/\" + merged_inhibitor_wins_name, format=\"sbml\")\n",
    "write_network_to_file(merged_network_or, home + \"LogicModelMerger/Models/\" + merged_or_name, format=\"sbml\")\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}