#!/usr/bin/env python
"""Provides functions to investigate the model and test with MEMOTE
These functions enable simple testing of any model using MEMOTE and access to its number of reactions, metabolites and genes.
"""
__author__ = "Famke Baeuerle and Alina Renz and Carolin Brune"
################################################################################
# requirements
################################################################################
import cobra
from cobra import Model as cobraModel
from libsbml import Model as libModel
from memote.support import consistency
# needed by memote.support.consistency
from memote.support import consistency_helpers as con_helpers
from ..developement.optional import require_bioregistry
from ..utility.util import test_biomass_presence
################################################################################
# variables
################################################################################
################################################################################
# functions
################################################################################
# get basic model info
# --------------------
[docs]
def get_reac_with_gpr(model: cobra.Model) -> tuple[list[str], list[str]]:
"""Extract the reactions that have a gene production rule from a given model.
Args:
- model (cobra.Model):
The model loaded with COBRApy.
Returns:
tuple:
Lists of reactions with GPR (1) & (2):
(1) list: List of normal reactions with gpr
(2) list: List of pseudoreactions with gpr
"""
pseudoreactions = model.boundary
for biomass in test_biomass_presence(model):
pseudoreactions.append(model.reactions.get_by_id(biomass))
normal = []
pseudo = []
for reac in model.reactions:
# check for GPR
if len(reac.genes) > 0:
if reac in pseudoreactions:
pseudo.append(reac.id)
else:
normal.append(reac.id)
return (normal, pseudo)
[docs]
def get_orphans_deadends_disconnected(
model: cobraModel,
) -> tuple[list[str], list[str], list[str]]:
"""Uses MEMOTE functions to extract orphans, deadends and disconnected metabolites
Args:
- model (cobraModel):
Model loaded with COBRApy
Returns:
tuple:
Lists of metabolites that might cause errors (1) - (3)
(1) list: List of orphans
(2) list: List of deadends
(3) list: List of disconnected metabolites
"""
orphans = consistency.find_orphans(model)
deadends = consistency.find_deadends(model)
disconnected = consistency.find_disconnected(model)
orphan_list = []
if len(orphans) > 0:
for orphan in orphans:
orphan_list.append(orphan.id)
deadend_list = []
if len(deadends) > 0:
for deadend in deadends:
deadend_list.append(deadend.id)
disconnected_list = []
if len(disconnected) > 0:
for disc in disconnected:
disconnected_list.append(disc.id)
return orphan_list, deadend_list, disconnected_list
[docs]
def get_mass_charge_unbalanced(model: cobraModel) -> tuple[list[str], list[str]]:
"""Creates lists of mass and charge unbalanced reactions, without exchange reactions since they are unbalanced per definition
Args:
- model (cobraModel):
Model loaded with COBRApy
Returns:
tuple:
Lists of reactions that might cause errors (1) & (2)
(1) list: List of mass unbalanced reactions
(2) list: List of charge unbalanced reactions
"""
mass_unbalanced = consistency.find_mass_unbalanced_reactions(model.reactions)
charge_unbalanced = consistency.find_charge_unbalanced_reactions(model.reactions)
mass_list = []
if len(mass_unbalanced) > 0:
for reac in mass_unbalanced:
if reac.id not in model.exchanges:
mass_list.append(reac.id)
charge_list = []
if len(charge_unbalanced) > 0:
for reac in charge_unbalanced:
if reac.id not in model.exchanges:
charge_list.append(reac.id)
return mass_list, charge_list
[docs]
def summarise_annotation_coverage(model:cobra.Model) -> tuple[dict, dict, dict]:
"""Summarise the annotation coverage of model entities by collecting types of annotations and their counts and percentages.
Args:
- model (cobra.Model): A model loaded with the COBRApy package.
Returns:
tuple[dict, dict, dict]:
Three dictionaries containing the annotation coverage for genes, reactions, and metabolites.
Each dictionary maps a database name to a tuple (count, percentage), where:
- database_name (str): name of the annotation database
- count (int): number of entities annotated with that database
- percentage (float): percentage of entities annotated with that database
.. code-block::
{
'database_name': (count, percentage),
...
}
"""
# genes
annot_genes = dict()
for gene in model.genes:
for db, _ in gene.annotation.items():
if db not in annot_genes:
annot_genes[db] = 0
annot_genes[db] += 1
total_genes = len(model.genes)
annot_genes = {db: (count, count/total_genes*100) for db, count in annot_genes.items()}
# reactions
annot_reacs = dict()
for reac in model.reactions:
for db, _ in reac.annotation.items():
if db not in annot_reacs:
annot_reacs[db] = 0
annot_reacs[db] += 1
total_reacs = len(model.reactions)
annot_reacs = {db: (count, count/total_reacs*100) for db, count in annot_reacs.items()}
# metabolites
annot_metabs = dict()
for metab in model.metabolites:
for db, _ in metab.annotation.items():
if db not in annot_metabs:
annot_metabs[db] = 0
annot_metabs[db] += 1
total_metabs = len(model.metabolites)
annot_metabs = {db: (count, count/total_metabs*100) for db, count in annot_metabs.items()}
return annot_genes, annot_reacs, annot_metabs
# other
# -----
# SBO terms
# ---------
[docs]
def get_reactions_per_sbo(model: libModel) -> dict:
"""Counts number of reactions of all SBO Terms present
Args:
- model (libModel):
Model loaded with libSBML
Returns:
dict:
SBO Term as keys and number of reactions as values
"""
is_valid_curie = require_bioregistry("validate SBO CURIEs").is_valid_curie
sbos_dict = {}
for react in model.getListOfReactions():
sbo = react.getSBOTerm()
complete_sbo_term = f'sbo:0000{sbo}'
if not is_valid_curie(complete_sbo_term):
sbo = 'invalid'
if sbo in sbos_dict.keys():
sbos_dict[sbo] += 1
else:
sbos_dict[sbo] = 1
return sbos_dict