GECKO
GECKO algorithm can be used to easily adjust the metabolic activity within the cell to respect many known parameters, measured by proteomics and other methods.
The original description from GECKO is by: Sánchez, et. al., "Improving the phenotype predictions of a yeast genome‐scale metabolic model by incorporating enzymatic constraints.", Molecular systems biology, 2017.
The analysis method and implementation in COBREXA is similar to sMOMENT, but GECKO is able to process and represent much larger scale of the constraints – mainly, it supports multiple isozymes for each reaction, and the isozymes can be grouped into "enzyme mass groups" to simplify interpretation of data from proteomics.
For demonstration, we will generate artificial random data in a way similar to the sMOMENT example:
!isfile("e_coli_core.json") &&
download("http://bigg.ucsd.edu/static/models/e_coli_core.json", "e_coli_core.json")
using COBREXA, GLPK
model = load_model("e_coli_core.json")
import Random
Random.seed!(1) # repeatability
gene_product_masses = Dict(genes(model) .=> randn(n_genes(model)) .* 10 .+ 60)
rxns = filter(
x ->
!looks_like_biomass_reaction(x) &&
!looks_like_exchange_reaction(x) &&
!isnothing(reaction_gene_association(model, x)),
reactions(model),
)69-element Vector{String}:
"PFK"
"PFL"
"PGI"
"PGK"
"PGL"
"ACALD"
"AKGt2r"
"PGM"
"PIt2r"
"ALCD2x"
⋮
"MALt2_2"
"MDH"
"ME1"
"ME2"
"NADH16"
"NADTRHD"
"NH4t"
"O2t"
"PDH"The main difference from sMOMENT comes from allowing multiple isozymes per reaction (reactions with missing isozyme informations will be ignored, leaving them as-is):
rxn_isozymes = Dict(
rxn => [
Isozyme(
Dict(isozyme_genes .=> 1),
randn() * 100 + 600, #forward kcat
randn() * 100 + 500, #reverse kcat
) for isozyme_genes in reaction_gene_association(model, rxn)
] for rxn in rxns
)Dict{String, Vector{Isozyme}} with 69 entries:
"ACALD" => [Isozyme(Dict("b0351"=>1), 462.073, 330.922), Isozyme(Dict("b124…
"PTAr" => [Isozyme(Dict("b2297"=>1), 471.533, 631.132), Isozyme(Dict("b245…
"ALCD2x" => [Isozyme(Dict("b0356"=>1), 617.367, 537.297), Isozyme(Dict("b124…
"PDH" => [Isozyme(Dict("b0114"=>1, "b0115"=>1, "b0116"=>1), 670.35, 423.7…
"PYK" => [Isozyme(Dict("b1676"=>1), 607.735, 454.975), Isozyme(Dict("b185…
"CO2t" => [Isozyme(Dict("s0001"=>1), 847.517, 338.625)]
"MALt2_2" => [Isozyme(Dict("b3528"=>1), 383.908, 397.408)]
"CS" => [Isozyme(Dict("b0720"=>1), 750.383, 549.73)]
"PGM" => [Isozyme(Dict("b0755"=>1), 586.875, 439.579), Isozyme(Dict("b361…
"TKT1" => [Isozyme(Dict("b2465"=>1), 650.793, 528.251), Isozyme(Dict("b293…
"ACONTa" => [Isozyme(Dict("b0118"=>1), 621.458, 449.334), Isozyme(Dict("b127…
"GLNS" => [Isozyme(Dict("b1297"=>1), 411.219, 513.669), Isozyme(Dict("b387…
"ICL" => [Isozyme(Dict("b4015"=>1), 657.74, 627.894)]
"FBA" => [Isozyme(Dict("b1773"=>1), 575.634, 511.58), Isozyme(Dict("b2097…
"FORt2" => [Isozyme(Dict("b0904"=>1), 669.647, 627.895), Isozyme(Dict("b249…
"G6PDH2r" => [Isozyme(Dict("b1852"=>1), 737.135, 408.509)]
"AKGDH" => [Isozyme(Dict("b0727"=>1, "b0116"=>1, "b0726"=>1), 537.049, 670.…
"TKT2" => [Isozyme(Dict("b2465"=>1), 676.15, 452.057), Isozyme(Dict("b2935…
"FRD7" => [Isozyme(Dict("b4151"=>1, "b4153"=>1, "b4154"=>1, "b4152"=>1), 5…
⋮ => ⋮We also construct similar bounds for total gene product amounts:
gene_product_bounds = Dict(genes(model) .=> Ref((0.0, 10.0)))Dict{String, Tuple{Float64, Float64}} with 137 entries:
"b4301" => (0.0, 10.0)
"b1602" => (0.0, 10.0)
"b4154" => (0.0, 10.0)
"b3236" => (0.0, 10.0)
"b1621" => (0.0, 10.0)
"b1779" => (0.0, 10.0)
"b3951" => (0.0, 10.0)
"b1676" => (0.0, 10.0)
"b3114" => (0.0, 10.0)
"b1241" => (0.0, 10.0)
"b2276" => (0.0, 10.0)
"b1761" => (0.0, 10.0)
"b3925" => (0.0, 10.0)
"b3493" => (0.0, 10.0)
"b3733" => (0.0, 10.0)
"b2926" => (0.0, 10.0)
"b0979" => (0.0, 10.0)
"b4015" => (0.0, 10.0)
"b2296" => (0.0, 10.0)
⋮ => ⋮With this, the construction of the model constrained by all enzymatic information is straightforward:
gecko_model =
model |> with_gecko(;
reaction_isozymes = rxn_isozymes,
gene_product_bounds,
gene_product_molar_mass = gene_product_masses,
gene_product_mass_group = _ -> "uncategorized", # all products belong to the same "uncategorized" category
gene_product_mass_group_bound = _ -> 100.0, # the total limit of mass in the single category
)Metabolic model of type GeckoModel sparse([9, 51, 55, 64, 65, 73, 9, 51, 55, 64 … 200, 201, 202, 203, 204, 205, 206, 207, 208, 209], [1, 1, 1, 1, 1, 1, 2, 2, 2, 2 … 325, 326, 327, 328, 329, 330, 331, 332, 333, 334], [1.0, 1.0, -1.0, -1.0, 1.0, -0.001546133025849219, 1.0, 1.0, -1.0, -1.0 … 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], 209, 334) Number of reactions: 334 Number of metabolites: 209
(Alternatively, you may use make_gecko_model, which does the same without piping by |>.)
The stoichiometry and coupling in the gecko model is noticeably more complex; you may notice new "reactions" added that simulate the gene product utilization:
[stoichiometry(gecko_model); coupling(gecko_model)]262×334 SparseArrays.SparseMatrixCSC{Float64, Int64} with 1386 stored entries:
⠄⠦⠤⠤⠄⠶⠶⠦⠇⠶⡴⠠⠀⠀⢧⠀⠢⣲⠼⠿⠶⣶⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⡎⢶⡉⣁⡄⣬⣠⡅⠔⠰⡈⠀⠀⠐⡄⡀⠄⣫⡉⡉⠩⠝⠇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠭⠰⠀⠒⣆⣀⣊⣬⢗⣾⢠⣵⣆⠆⠀⡁⠨⠁⠈⠀⠂⠥⠁⠅⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⣧⡣⠁⠤⠟⠒⠻⠿⡗⠆⠠⣧⣤⢵⠀⢠⣄⡀⠘⠋⠚⠠⠂⠆⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⢶⣂⠀⠀⠀⠀⠀⠀⠃⠀⠀⠉⠉⠉⠁⠈⠉⠉⠓⠀⠀⠀⠀⠰⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠈⠙⠺⠤⣀⡀⠀⠤⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠤⠀⠠⠄⠀⠱⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠘⣆⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠁⠀⠀⠀⠑⢆⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠉⠳⢄⡀⠀⠀⠀⠀⠀⠐⠀⠀⠀⠐⠀⠀⠀⠀⠀⠀⠀⠑⢆⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⠲⠄⠀⠠⣄⠀⠀⠀⠀⠀⠐⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠧⣜⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀
⠒⠰⠤⣄⣀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⠃⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙
⠀⠀⠀⠀⠀⠉⠁⠲⠤⢤⡀⣀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠙⠀⠀⠐⠲⠤⢠⣀⣀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⠘⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒Again, the resulting model can be used in any type of analysis. For example, flux balance analysis:
opt_model = flux_balance_analysis(gecko_model, GLPK.Optimizer)A JuMP Model
Maximization problem with:
Variables: 334
Objective function type: JuMP.AffExpr
`JuMP.AffExpr`-in-`MathOptInterface.EqualTo{Float64}`: 209 constraints
`JuMP.AffExpr`-in-`MathOptInterface.LessThan{Float64}`: 774 constraints
Model mode: AUTOMATIC
CachingOptimizer state: ATTACHED_OPTIMIZER
Solver name: GLPK
Names registered in the model: c_lbs, c_ubs, lbs, mb, ubs, xGet the fluxes
flux_sol = flux_dict(gecko_model, opt_model)Dict{String, Float64} with 95 entries:
"ACALD" => 0.0
"PTAr" => 6.86778
"ALCD2x" => 0.0
"PDH" => 10.5594
"PYK" => 2.58923
"CO2t" => -15.7718
"EX_nh4_e" => -3.90714
"MALt2_2" => 0.0
"CS" => 1.0062
"PGM" => -15.0145
"TKT1" => 1.8806
"EX_mal__L_e" => 0.0
"ACONTa" => 1.0062
"EX_pi_e" => -2.63593
"GLNS" => 0.183219
"ICL" => 0.0
"EX_o2_e" => -14.9434
"FBA" => 7.27848
"EX_gln__L_e" => 0.0
⋮ => ⋮Get the gene product concentrations
gp_concs = gene_product_dict(gecko_model, opt_model)Dict{String, Float64} with 137 entries:
"b4301" => 0.0
"b1602" => 0.0
"b4154" => 0.0
"b3236" => 0.000313395
"b1621" => 0.0
"b3951" => -7.43092e-19
"b1779" => 0.025239
"b1676" => -2.92291e-18
"b3114" => 0.0
"b1241" => 0.0
"b2276" => 0.0488441
"b1761" => 0.00556268
"b3925" => 0.0
"b3493" => 0.00343844
"b3733" => 0.0534869
"b2926" => 0.0332595
"b0979" => 0.0
"b4015" => 0.0
"b2296" => 0.0
⋮ => ⋮Get the total masses assigned to each mass group
gene_product_mass_group_dict(gecko_model, opt_model)Dict{String, Float64} with 1 entry:
"uncategorized" => 100.0Variability:
flux_variability_analysis(gecko_model, GLPK.Optimizer, bounds = gamma_bounds(0.95))95×2 Matrix{Float64}:
6.05514 11.1416
0.0 14.0378
0.370962 9.86045
-18.1905 -14.2797
1.61973e-15 9.48949
-3.37813 0.0
-1.16712 0.0
-17.1721 -13.2614
2.50413 2.63593
-3.37813 0.0
⋮
0.0 0.0
-3.14082 7.2849
0.0 2.97194
0.0 4.02313
26.5111 30.5831
0.0 7.30744
3.71178 4.72073
13.7438 16.3803
0.0 15.4793...and sampling:
affine_hit_and_run(gecko_model, warmup_from_variability(gecko_model, GLPK.Optimizer))' *
reaction_flux(gecko_model)3340×95 Matrix{Float64}:
5.69741 7.14345 4.35447 -9.7528 … 0.179103 4.165 1.85149
6.17793 7.68231 4.74709 -9.98991 0.161595 3.95987 1.48158
6.54096 8.12222 4.88527 -10.4291 0.130186 4.47108 1.55276
5.79829 7.23795 4.41735 -9.92733 0.179442 4.31659 1.8863
5.70484 7.46217 4.27064 -9.51379 0.118997 4.03253 1.63139
5.7485 7.1552 4.41012 -9.83999 … 0.187914 4.22067 1.89546
5.72688 7.18375 4.36675 -9.76999 0.174085 4.17519 1.83555
5.52281 6.53495 4.13366 -9.68977 0.251199 4.17018 2.46207
5.71168 7.20184 4.38244 -9.73635 0.172702 4.04388 1.8206
5.6307 7.0402 4.34134 -9.58676 0.154022 3.96686 1.89444
⋮ ⋱
5.7591 7.3333 4.18016 -9.80307 0.0971256 3.92446 2.01831
6.11556 7.58397 4.89516 -10.3363 0.186786 3.83085 1.8704
6.62928 7.69967 5.02136 -10.9306 0.178771 4.37013 2.54907
6.0239 7.6031 4.70331 -10.3331 0.172041 4.01933 2.01269
5.99965 7.59361 4.683 -10.3152 … 0.171666 4.04386 2.01911
5.92698 7.97307 4.70791 -10.2445 0.199597 4.04475 1.55914
5.97615 8.07602 4.67444 -10.3848 0.154037 4.17861 1.77528
6.13383 7.96383 5.18158 -10.8211 0.175543 3.28189 1.88694
5.65149 7.75653 4.44149 -9.97254 0.150216 4.02244 1.70548This page was generated using Literate.jl.