itle: Difficulty Performing SteadyCom Analysis on Genome-Scale Community Model Created with CarveMe

[wtscott31]

Hello @cdanielmachado, other developers and all in the community,

I am interested in performing analysis on my genome-scale community model, which I created using CarveMe. Specifically, I want to utilize some of the functions found in ReFramed, including the SteadyCom function. However, I am facing two main issues:

Lack of Documentation: I could not find much documentation on how to use the functions provided by ReFramed, including SteadyCom. Detailed examples and explanations would be extremely helpful.

KeyError: I keep encountering the following error when I attempt to run my analysis:

KeyError: 'community_biomass_dehalobacter_community_carveme'

Here is the code I am using:

import reframed
from reframed import load_cbmodel
from reframed.community.model import Community
from reframed.community.SteadyCom import SteadyCom
import pandas as pd
import matplotlib.pyplot as plt

# Path to the community model SBML file
community_model_path = "~/dehalobacter_community_carveme.xml"

# Load the community model
community_model = load_cbmodel(community_model_path)

# Initialize Community object with the loaded community model
community = Community("dehalobacter_community", [community_model])

# Define a new medium for the community in ReFramed format
new_medium_reframed = {
    'EX_meoh_e': 10.0,
    'EX_124triclbenz_e': 10.0,
    'EX_co2_e': 50.0,
    'EX_ca2_e': 1000,
    'EX_fe2_e': 1000,
    'EX_mobd_e': 1000,
    'EX_thm_e': 0.1,
    'EX_fol_e': 0.1,
    'EX_cl_e': 1000.0,
    'EX_fe3_e': 0.1,
    'EX_h2_e': 1000.0,
    'EX_h_e': 10,
    'EX_h2o_e': 1000.0,
    'EX_k_e': 1000.0,
    'EX_mg2_e': 1000.0,
    'EX_nh4_e': 1000,
    'EX_pi_e': 1000.0,
    'EX_so4_e': 1000.0,
    'EX_mn2_e': 1000.0,
    'EX_zn2_e': 1000.0,
    'EX_cobalt2_e': 1000.0,
    'EX_cu2_e': 1000.0,
    'EX_h2s_e': 5.0,
}

# Update the community's medium
community.medium = new_medium_reframed

# Perform SteadyCom analysis
steadycom_result = SteadyCom(community)

# Extract and print results
print(steadycom_result)

# Convert SteadyCom results to DataFrame
df_steadycom = pd.DataFrame.from_dict(steadycom_result.values, orient='index', columns=['flux'])
df_steadycom.reset_index(inplace=True)
df_steadycom.columns = ['reaction', 'flux']

# Save SteadyCom results to CSV
df_steadycom.to_csv('dehalobacter_community_steadycom_results.csv', index=False)

# Plot SteadyCom results
plt.figure(figsize=(10, 6))
df_steadycom_sorted = df_steadycom.sort_values(by='flux', ascending=False).head(20)  # Select top 20 reactions by flux
df_steadycom_sorted.plot(x='reaction', y='flux', kind='bar', legend=False, color='green')
plt.title('SteadyCom Analysis Results')
plt.ylabel('Flux')
plt.xlabel('Reaction')
plt.tight_layout()
plt.xticks(rotation=90)
plt.savefig('dehalobacter_community_steadycom_results.png')
plt.show()

It is odd that I get this KeyError since the model I have contains a biomass equation called "community_growth".

Could you please provide guidance on how to resolve the KeyError and any additional documentation or examples for using ReFramed's functions, particularly SteadyCom? I can also provide the genome-scale community model (GEM) file if needed.

Thank you for your assistance!

[cdanielmachado]

Hi, you are totally correct, I apologize for the lack of documentation... A lot of the functionality in this package is mostly experimental :)

CarveMe can generate an integrated community model as a single SBML file. But that model is only meant for a classical FBA simulation.

To take advantage of more sophisticated community simulation methods like SteadyCom, you need to provide a list of individual species models when you build the community.

That means that this:

# Initialize Community object with the loaded community model
community = Community("dehalobacter_community", [community_model])

Should be replaced with:

community = Community("dehalobacter_community", [species1, species2, ... speciesN])

I hope this helps :)

[wtscott31]

Hi @cdanielmachado

I revised my script to call each model separately for the SteadyCom function. Despite these adjustments, I am still encountering errors related to the Gurobi solver. Below are the details of the error and my current script implementation.

Error Details:

`AttributeError Traceback (most recent call last)
Input In [11], in <cell line: 60>()
57 type(community.solver)
59 # Perform SteadyCom analysis
---> 60 steadycom_result = SteadyCom(community)
61 # Extract and print results
62 print(steadycom_result)

File /opt/anaconda3/lib/python3.9/site-packages/reframed/community/SteadyCom.py:26, in SteadyCom(community, constraints, solver)
22 solver = build_problem(community)
24 objective = {community.merged_model.biomass_reaction: 1}
---> 26 sol = binary_search(solver, objective, minimize=False, constraints=constraints)
28 solution = CommunitySolution(community, sol.values)
29 solution.solver = solver

File /opt/anaconda3/lib/python3.9/site-packages/reframed/community/SteadyCom.py:156, in binary_search(solver, objective, obj_frac, minimize, max_iters, abs_tol, constraints)
153 fold = 0.5
154 value = fold*diff + previous_value
--> 156 solver.update_growth(value)
157 sol = solver.solve(objective, get_values=False, minimize=minimize, constraints=constraints)
159 feasible = sol.status == Status.OPTIMAL

File /opt/anaconda3/lib/python3.9/site-packages/reframed/community/SteadyCom.py:126, in build_problem..update_growth(value)
123 def update_growth(value):
124 # TODO: find a solution that is not CPLEX specific
125 coefficients = [(f"g_{x}", f"x_{x}", value) for x in community.organisms]
--> 126 solver.problem.linear_constraints.set_coefficients(coefficients)

File src/gurobipy/model.pxi:357, in gurobipy.Model.getattr()

File src/gurobipy/model.pxi:1901, in gurobipy.Model.getAttr()

File src/gurobipy/attrutil.pxi:23, in gurobipy.__getattrinfo()

AttributeError: 'gurobipy.Model' object has no attribute 'linear_constraints'
`

For the revised script:

import reframed
from reframed import load_cbmodel
from reframed.community.model import Community
from reframed.community.SteadyCom import SteadyCom
import pandas as pd
import matplotlib.pyplot as plt

# Path to the community model SBML file
Dehalobacter_restrictus = load_cbmodel("Gap-filled_MM2014_GEMs/Dehalobacter_restrictus_new_3.xml")
Dehalobacter_sp004343605 = load_cbmodel("Gap-filled_MM2014_GEMs/Dehalobacter_sp004343605_new_3.xml")
Acetobacterium_sp003260995 = load_cbmodel("Gap-filled_MM2014_GEMs/Acetobacterium_sp003260995.xml")
Cloacimonadaceae_g__Cloacimonas = load_cbmodel("Gap-filled_MM2014_GEMs/Cloacimonadaceae_g__Cloacimonas.xml")
g__Methanothrix_s = load_cbmodel("Gap-filled_MM2014_GEMs/g__Methanothrix_s.xml")
Methanobacterium_C_congolense = load_cbmodel("Gap-filled_MM2014_GEMs/Methanobacterium_C_congolense.xml")
o__Bacteroidales_f__TTA_H9_g__TTA_H9_s__TTA_H9_sp002418865 = load_cbmodel("Gap-filled_MM2014_GEMs/o__Bacteroidales_f__TTA-H9_g__TTA-H9_s__TTA-H9_sp002418865.xml")
o__Bacteroidales_f__VadinHA17_g__LD21_s__ = load_cbmodel("Gap-filled_MM2014_GEMs/o__Bacteroidales_f__VadinHA17_g__LD21_s__.xml")
o__Bacteroidales_f__VadinHA17_g__SR_FBR_E99_s__ = load_cbmodel("Gap-filled_MM2014_GEMs/o__Bacteroidales_f__VadinHA17_g__SR-FBR_E99_s__.xml")
Rectinema_sp002441395 = load_cbmodel("Gap-filled_MM2014_GEMs/Rectinema_sp002441395.xml")
Rectinema_sp015657395 = load_cbmodel("Gap-filled_MM2014_GEMs/Rectinema_sp015657395.xml")
Rectinema_subterraneum = load_cbmodel("Gap-filled_MM2014_GEMs/Rectinema_subterraneum.xml")

# Initialize Community object with the loaded community model
community = Community("dehalobacter_community", [Dehalobacter_restrictus, Dehalobacter_sp004343605, Acetobacterium_sp003260995, Cloacimonadaceae_g__Cloacimonas, g__Methanothrix_s, Methanobacterium_C_congolense, o__Bacteroidales_f__TTA_H9_g__TTA_H9_s__TTA_H9_sp002418865, o__Bacteroidales_f__VadinHA17_g__LD21_s__, o__Bacteroidales_f__VadinHA17_g__SR_FBR_E99_s__, Rectinema_sp002441395, Rectinema_sp015657395, Rectinema_subterraneum])

# Define a new medium for the community in ReFramed format
new_medium_reframed = {
    'EX_meoh_e': 10.0,
    'EX_124triclbenz_e': 10.0,
    'EX_co2_e': 50.0,
    'EX_ca2_e': 1000,
    'EX_fe2_e': 1000,
    'EX_mobd_e': 1000,
    'EX_thm_e': 0.1,
    'EX_fol_e': 0.1,
    'EX_cl_e': 1000.0,
    'EX_fe3_e': 0.1,
    'EX_h2_e': 1000.0,
    'EX_h_e': 10,
    'EX_h2o_e': 1000.0,
    'EX_k_e': 1000.0,
    'EX_mg2_e': 1000.0,
    'EX_nh4_e': 1000,
    'EX_pi_e': 1000.0,
    'EX_so4_e': 1000.0,
    'EX_mn2_e': 1000.0,
    'EX_zn2_e': 1000.0,
    'EX_cobalt2_e': 1000.0,
    'EX_cu2_e': 1000.0,
    'EX_h2s_e': 5.0,
}

# Update the community's medium
community.medium = new_medium_reframed

community.solver = 'glpk'
type(community.solver)

# Perform SteadyCom analysis
steadycom_result = SteadyCom(community)
# Extract and print results
print(steadycom_result)

# Convert SteadyCom results to DataFrame
df_steadycom = pd.DataFrame.from_dict(steadycom_result.values, orient='index', columns=['flux'])
df_steadycom.reset_index(inplace=True)
df_steadycom.columns = ['reaction', 'flux']

# Save SteadyCom results to CSV
df_steadycom.to_csv('dehalobacter_community_steadycom_results.csv', index=False)

# Plot SteadyCom results
plt.figure(figsize=(10, 6))
df_steadycom_sorted = df_steadycom.sort_values(by='flux', ascending=False).head(20)  # Select top 20 reactions by flux
df_steadycom_sorted.plot(x='reaction', y='flux', kind='bar', legend=False, color='green')
plt.title('SteadyCom Analysis Results')
plt.ylabel('Flux')
plt.xlabel('Reaction')
plt.tight_layout()
plt.xticks(rotation=90)
plt.savefig('dehalobacter_community_steadycom_results.png')
plt.show()

The SteadyCom analysis should run without encountering the AttributeError.

Additional Context

Operating System: macOS
Python Version: 3.9
ReFramed Version: 1.2.1
Gurobi Version: 1.5.1

Any guidance or suggestions to resolve this issue would be greatly appreciated.

Thank you!

[cdanielmachado]

I suppose this is the problem... 😅

def update_growth(value):
        # TODO: find a solution that is not CPLEX specific
        coefficients = [(f"g_{x}", f"x_{x}", value) for x in community.organisms]
        solver.problem.linear_constraints.set_coefficients(coefficients)

[wtscott31]

Hi @cdanielmachado Thanks for your reply! Yes, indeed the SteadyCom function works fine when using CPLEX as the solver.