add visualization of SBO-Terms #56

main.py

@@ -1,9 +1,11 @@
 #!/usr/bin/env python
 
 import yaml
+import os
 import refinegems as rg
 import cobra
 import pandas as pd
+import matplotlib.pyplot as plt
 from datetime import date
 
 __author__ = "Famke Baeuerle and Gwendolyn O. Gusak"
@@ -18,6 +20,13 @@ def main():
 
     with open('config.yaml') as f:
         config = yaml.safe_load(f)
+
+    # check if the output directory is already present, if not create it
+    dir = os.path.join(config['out_path'] + 'visualization/')
+    if not os.path.isdir(config['out_path']):
+        os.makedirs(config['out_path'])
+    if not os.path.isdir(dir): 
+        os.makedirs(dir)
 
     if (config['keggpathways']):
         non_kegg = rg.pathways.kegg_pathways(config['model'], config['kegg_path'])
@@ -63,6 +72,8 @@ def main():
         if (config['multiple']):
             growth_all = rg.comparison.simulate_all(config['multiple_paths'], config['media'], config['growth_basis'])
             growth_all.to_csv(config['out_path'] + 'growth_' + str(today) + '_' + config['growth_basis'] + '.csv', index=False)
+            sbo_fig_all = rg.comparison.get_sbo_plot_multiple(config['multiple_paths']).get_figure()
+            sbo_fig_all.savefig(config['out_path'] + 'visualization/' + 'all_ReacPerSBO_' + str(today) + '.png', bbox_inches='tight')
 
         try:    
             model_cobra, errors = cobra.io.sbml.validate_sbml_model(config['model'])
@@ -77,6 +88,10 @@ def main():
             orphans, deadends, disconnected = rg.investigate.get_orphans_deadends_disconnected(model_cobra)
             mass_unbal, charge_unbal = rg.investigate.get_mass_charge_unbalanced(model_cobra)
             egc = rg.investigate.get_egc(model_cobra)
+            sbo_fig = rg.investigate.get_sbo_plot_single(model_libsbml).get_figure()
+
+            # saving the created visualizations
+            sbo_fig.savefig(config['out_path'] + 'visualization/' + str(model_cobra.id) + '_ReacPerSBO_' + str(today) + '.png', bbox_inches='tight')
 
             if (config['memote']):
                 score = rg.investigate.run_memote(model_cobra)

refinegems/comparison.py

@@ -6,40 +6,38 @@
 
 import pandas as pd
 import matplotlib.pyplot as plt
-from cobra import Model
+from cobra import Model as cobmod
+from libsbml import Model as libmod
 from tqdm import tqdm
 from venn import venn
-from refinegems.io import load_multiple_models, load_all_media_from_db
+from refinegems.io import load_multiple_models, load_medium_from_db, search_sbo_label
 from refinegems.growth import growth_one_medium_from_default, growth_one_medium_from_minimal
-from refinegems.investigate import initial_analysis
+from refinegems.investigate import initial_analysis, get_reactions_per_sbo
 
 __author__ = "Famke Baeuerle"
 
-sbo_mapping={'658': 'passive transport', 
-            '176': 'biochemical reaction', 
-            '167': 'biochemical or transport reaction',
-            '402': 'transfer of a chemical group', 
-            '659': 'symporter-mediated transport', 
-            '200': 'redox reaction', 
-            '233': 'hydroxylation',
-            '399': 'decarboxylation', 
-            '178': 'cleavage', 
-            '403': 'transamination', 
-            '215': 'acetylation', 
-            '377': 'isomerisation', 
-            '657': 'active transport', 
-            '216': 'phosphorylation', 
-            '401': 'deamination', 
-            '376': 'hydrolysis', 
-            '217': 'glycosylation', 
-            '660': 'antiporter-mediated transport', 
-            '654': 'co-transport reaction', 
-            '214': 'methylation', 
-            '655': 'transport reaction', 
-            '627': 'exchange reaction', 
-            '632': 'sink reaction', 
-            '629': 'biomass production',
-            '630': 'ATP maintenance'}
+def get_sbo_mapping_multiple(models):
+    mappings = {}
+    for model in models:
+        mappings[model.id] = get_reactions_per_sbo(model)
+    df = pd.DataFrame.from_dict(mappings)
+    df = df.reset_index().rename({'index': 'SBO-Term'}, axis=1)
+    df['SBO-Name'] = df['SBO-Term'].apply(search_sbo_label)
+    return df
+
+def get_sbo_plot_multiple(model_list: list[str], rename=None):
+    models = load_multiple_models(model_list, package='libsbml')
+    map = get_sbo_mapping_multiple(models)
+    id_list = [mod.id for mod in models]
+    map = map[(map[id_list]>3).all(axis=1)]
+    map = map.drop('SBO-Term', axis=1).sort_values(id_list[0]).set_index('SBO-Name')
+    if rename is not None:
+        map = map.rename(rename, axis=1)
+    fig = map.plot.barh(stacked=True, width=.8, figsize=(8,10))
+    fig.set_ylabel('')
+    fig.set_xlabel('number of reactions', fontsize=16)
+    fig.legend(loc='lower right')
+    return fig
 
 def create_venn(model_list: list[str], entity: str, perc: bool=False) -> plt.figure: 
     all_models = load_multiple_models(model_list, package='cobra')
@@ -59,7 +57,7 @@ def create_venn(model_list: list[str], entity: str, perc: bool=False) -> plt.fig
         fig = venn(intersec)
     return fig
 
-def simulate_all(model_list: list[str], mediumpath: str, media: list[str], basis: str) -> pd.DataFrame:
+def simulate_all(model_list: list[str], media: list[str], basis: str) -> pd.DataFrame:
     """does a run of growth simulation for multiple models on different media
 
     Args:
@@ -72,10 +70,9 @@ def simulate_all(model_list: list[str], mediumpath: str, media: list[str], basis
         df: table containing the results of the growth simulation
     """
     growth = pd.DataFrame()
-    all_media = load_all_media_from_db(mediumpath)
     all_models = load_multiple_models(model_list, package='cobra')
-    selected_media = [x for x in all_media if x['medium'][0] in media]
-    for medium in tqdm(selected_media):
+    for medium_id in tqdm(media):
+        medium = load_medium_from_db(medium_id)
         for model in all_models:
             essentials_given = False
             if (basis=='default_uptake'):

refinegems/investigate.py

@@ -14,7 +14,7 @@
 from memote.support import consistency
 # needed by memote.support.consitency
 from memote.support import consistency_helpers as con_helpers
-from refinegems.io import load_model_cobra, load_model_libsbml
+from refinegems.io import load_model_cobra, load_model_libsbml, search_sbo_label
 
 __author__ = "Famke Baeuerle"
 
@@ -283,4 +283,24 @@ def get_metabs_with_one_cvterm(model):
             if sb.getCVTerm(0).getNumResources() == 1:
                 only_one.append(pid)
 
-    return only_one
+    return only_one
+
+def get_reactions_per_sbo(model):
+    sbos_dict = {}
+    for react in model.getListOfReactions():
+        sbo = react.getSBOTerm()
+        if sbo in sbos_dict.keys():
+            sbos_dict[sbo] += 1
+        else: 
+            sbos_dict[sbo] = 1
+    return sbos_dict
+
+def get_sbo_plot_single(model):
+    df = pd.DataFrame(get_reactions_per_sbo(model), index=[0]).T.reset_index().rename({0:model.id, 'index': 'SBO-Term'}, axis=1)
+    df = df[df[model.id]>3]
+    df['SBO-Name'] = df['SBO-Term'].apply(search_sbo_label)
+    fig = df.drop('SBO-Term', axis=1).sort_values(model.id).set_index('SBO-Name').plot.barh(width=.8, figsize=(8,10))
+    fig.set_ylabel('')
+    fig.set_xlabel('number of reactions', fontsize=16)
+    fig.legend(loc='lower right')
+    return fig

refinegems/io.py

@@ -11,6 +11,7 @@
 import gffutils
 import sqlalchemy
 import pandas as pd
+from ols_client import EBIClient
 from Bio import Entrez, SeqIO
 from refinegems.databases import PATH_TO_DB
 from libsbml import SBMLReader, writeSBMLToFile, Model, SBMLValidator, SBMLDocument
@@ -328,4 +329,10 @@ def extract_locus(feature):
 
     mapping_df['locus_tag'] = mapping_df.apply(
         lambda row: extract_locus(row['Parent']), axis=1)
-    return mapping_df.drop('Parent', axis=1)
+    return mapping_df.drop('Parent', axis=1)
+
+def search_sbo_label(sbo_number: str):
+    sbo_number = str(sbo_number)
+    client = EBIClient()
+    sbo = client.get_term('sbo', 'http://biomodels.net/SBO/SBO_0000' + sbo_number)
+    return sbo['_embedded']['terms'][0]['label']